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SUMMARY 

A detailed  description  is  given  of  the  computer  program  SYNMAP,  which  uses 
a stochastic  simulation  method  to  determine  the  total  velocity  increment 
required  for  the  station  acquisition  phase  of  a synchronous  satellite  orbit 
mission.  The  program  takes  account  of  errors  due  to  launch  vehicle  injection, 
satellite  tracking  and  apogee  motor  burn.  A description  is  also  given  of  the 
program  P0INT2,  which  may  be  used  to  generate  the  set  of  random  transfer  orbits 
required  by  SYNMAP. 
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1 INTRODUCTION 

The  program  SYNMAP  (an  acronym  for  synchronous  mission  analysis  £rogram) 

is  part  of  the  computer  software  required  for  analysing  the  launch  phase  of  a 

synchronous  satellite  mission.  The  program  uses  a stochastic  simulation  method 

to  determine  the  total  velocity  increment  required  for  station  acquisition, 

taking  account  of  launch  vehicle  injection  errors  and  apogee  motor  burn  errors. 

Tracking  errors  may  also  be  included,  but  their  effect  is  small  and  they  are 

normally  neglected.  The  program  may  be  used  for  any  synchronous  mission  in 

which  the  spacecraft  is  injected  into  a transfer  orbit  and  then  inserted  into  a 

drift  orbit  by  use  of  a fixed  impulse  apogee  motor.  Missions  for  which  the 

1 2 

program  has  been  used  include  SKYNET  3 and  MAROTS  . 

The  program  must  be  supplied  with  a set  of  random  transfer  orbits,  each  of 
which  is  specified  by  the  cartesian  components  of  the  geocentric  position  and 
velocity  vectors  at  a time  shortly  before  the  nominal  apogee  boost  motor  (ABM) 
firing  point.  These  orbits  are  usually  generated  by  the  program  P0INT2,  a modi- 
fied version  of  POINT  which  is  described  in  detail  in  Ref  3. 

SYNMAP  provides  a choice  of  strategy  for  both  ABM  firing  and  station 
acquisition.  The  firing  strategy  may  optimise  the  drift  orbit  flight  path  angle 
to  minimise  the  total  velocity  increment  required  for  station  acquisition, 
provide  a given  flight  path  angle  or  provide  a given  drift  rate.  The  station 
acquisition  strategy  enables  the  satellite  to  be  placed  on  station  with  either 
an  absolute  minimum  velocity  increment  or  a minimum  consistent  with  the  satel- 
lite being  moved  through  the  smaller  longitude  range. 

The  output  from  SYNMAP  consists  of  two  lines  for  each  simulation  and  a 
block  of  statistical  information.  For  each  simulation,  the  program  gives  the 
right  ascension  and  declination  of  the  ABM  thrust  direction,  the  drift  orbit 
osculating  elements  immediately  after  ABM  burn,  the  longitude  of  the  burn  point 
and  the  delta-velocity  components  required  for  station  acquisition  and  circular- 
ization of  the  orbit.  The  statistical  information  consists  of  the  mean  and 
standard  deviation  for  the  following  quantities:  drift  orbit  semi  major  axis, 
inclination,  eccentricity  and  right  ascension  of  the  ascending  mode;  ABM  firing 
longitude;  right  ascension  and  declination  of  the  ABM  thrust  direction;  solar 
aspect  angle  and  the  total  velocity  increment  required  for  station  acquisition. 

A table  giving  a histogram  of  the  velocity  increments  required  is  also  printed. 

The  software  is  written  in  1900  FORTRAN  for  use  on  ICL  1900  series  com- 
puters. P0INT2  requires  approximately  I5K  words  of  core  store,  and  consists  of 
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a main  program  and  18  aubprograms,  while  SYNMAP  requires  approximately  20K  words 
of  core  store,  and  consists  of  a main  program  and  32  subprograms.  The  program 
units  and  subprogram  specifications  for  P0INT2  are  given  in  Appendices  A and  B, 
and  those  for  SYNMAP  in  Appendices  C and  D.  Their  calling  structures  are 
illustrated  in  Figs  I and  2. 

2 DESCRIPTION  OF  P0IKT2 

The  program  P0INT2  is  used  to  generate  two  sets  of  random  transfer  orbits 
for  use  with  SYNMAP.  One  set  is  for  the  prime  ABM  firing  apogee  and  the  second 
for  the  back-up  apogee.  Each  set  has  a common  epoch,  usually  three  hours  before 
apogee,  and  each  orbit  is  defined  by  its  geocentric  position  and  velocity  vectors 
at  epoch. 

The  program  is  provided  with  a set  of  nominal  launch  vehicle  injection  con- 
ditions (velocity,  climb  angle,  azimuth,  radial  distance,  latitude  and  longi- 
tude). These,  represented  by  the  elements  of  the  column  vector  x , are  assumed 
to  have  a multivariate  normal  distribution  about  their  mean  vector  £ , such 
that  the  probability  density  function  of  x is  given  by 

f(x)  - £(2ir)n|v|J  exp|~-  i(x“£)TV  1 (x  " w)J  O) 

where  V is  the  covariance  error  matrix  and  n is  the  number  of  variates  (six 
in  the  present  case). 

Simulation  of  the  injection  process  requires  the  generation  of  random 

vectors  from  the  distribution  defined  by  equation  (1);  the  stochastically 

dependent  elements  of  this  vector  can  be  constructed  from  a set  of  independent 
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normal  variables  using  the  linear  transformation  method  suggested  by  Moonan  . 
Consider  a column  vector  £ of  m independent  variables  (z^,  i = |,...,m) 


E(z.)  - 0 , 


1 < i < m 


E(zz;  ■ 1^  , 

where  E is  expectation  and  1^  is  the  unit  matrix  of  order  m x m . _z  may 
be  transformed  using  the  linear  mapping 

x - £ ■ Cz  , (2) 

where  C is  a lower-triangular  nonsingular  matrix  of  order  m x m.  For  the  new 
vector  x Co  have  the  probability  density  function  of  equation  (1), 


E(x.) 


Vi 


1 < i < m 


E[(x  - y)(x  - y)TJ  - e|czzTCT1  - CE(zzT)CT  = CC1 


- V 


(3) 


Hence  to  determine  the  transformation,  the  symmetric  and  positive  definite 


matrix  V is  decomposed  into  the  product  of  a lower  triangular  matrix  C and 
. T 

its  transpose  C . The  elements  (c..)  of  C can  be  determined  recursively 


from  the  equation. 


ij 


Cil 


v. 


• /v  1 

il/  11 


1 £ i ^ m 


i-1 


c. . 
il 


r V 2 y 

■ V. . - ) C.,  , 

L ii  Lj  ikJ 


2 ^ i < m 


k-1 

i-1 


>►  (4) 


c. . 

iJ 


[Vij  ^ CikCjk]/Cjj  * 


1 < j < i ^ m 


k-1 


and 


c. . 
U 


- 0 , 


J > i 


J 


To  obtain  each  random  vector  x , six  independent  normal  variates,  with 
zero  mean  and  unit  variance,  are  generated  and  the  transformation  of  equation  (2) 


applied.  If  r^  and  r^  are  independent  random  variables  from  a uniform  dis- 


tribution defined  on  the  interval  (0,1),  then 


(-  2 In  Tj)5  cos  2^2 


and 


x2  ■ (-  2 In  sin  2irr2 


(5) 


are  a pair  of  independent  random  variables  from  a normal  distribution  with  zero 
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mean  and  unit  variance  . The  program  uses  equation  (5)  with  pseudo-random  num- 
bers for  Tj  and  r produced  by  the  mixed  congruent ial  method,  ie  using 


generators  of  the  form 


n^+j  - an^  + c (mod  m)  , 


(6) 


where  n. , a,  c and  m are  all  non-negative  integers. 


A set  of  geocentric  position  and  velocity  components  is  formed  from  each 
vector  x and  the  orbit  integrated  to  the  first  epoch  using  the  method 
described  in  section  3.  An  interpolation  is  performed  and  the  position  and 
velocity  components  are  written  to  a disc  file.  The  integration  is  then  con- 
tinued to  the  second  epoch  and  the  interpolated  position  and  velocity  components 
written  to  a second  disc  file. 

If  required,  the  integration  may  be  halted  at  a given  apogee  and  a gross 
attitude  manoeuvre  simulated  by  adding  the  appropriate  velocity  increments  to 
the  velocity  vector.  The  integration  is  then  restarted  and  continued  as  before. 

To  enable  sets  of  orbits  to  be  reproduced,  the  random  number  generator  is 
initialised  using  constants  stored  within  the  program. 

3 INTEGRATION  PROCEDURE 

The  equation  of  motion  of  an  earth  satellite  may  be  expressed  in  vector 
form  as 

-3 

£ - - pr  r + F (7) 

where  £ is  the  position  vector  relative  to  the  earth's  centre  of  mass, 
p (»  GM)  is  the  earth's  gravitational  constant  and  F is  the  perturbing 
acceleration.  Orbit  development  is  obtained  by  the  numerical  integration  of  the 
equations  of  motion  as  formulated  in  the  Cowell  method.  Because  of  the  large 
variations  in  the  velocity  and  perturbing  accelerations  around  a highly  eccentric 
orbit,  a constant  time  interval  cannot  be  used  if  the  computation  is  to  be 
efficient.  To  overcome  this  difficulty,  the  equations  are  transformed  in  such 
a way  that  a constant  step  length  may  be  used,  ie  analytical  step  regulation  is 
introduced.  Time  t is  replaced  by  the  independent  variable  s , defined  by 

dt/ds  ■ r^/x  (8) 

where  x is  a constant.  Merson  has  discussed^  the  selection  of  k and  x to 
give  the  best  results  and  suggested  the  use  of  k ■ 3/2  and  x = these 
values  are  used  here.  On  changing  to  the  independent  variable  s , where 


w*wiiiw}WBWWtiBP  W-  [iMipi!  ijiirmr.i-ijijfir 


The  equation  for  t'  can  be  differentiated  to  give 


t"  “ I (r  . r')/(ur)^ 


(ID 


so  that  we  have  four  second-order  differential  equations  (for  x",  y",  z",  t"). 
The  numerical  integration  in  both  P0INT2  and  SYNMAP  is  based  on  an  eighth-order 
Gauss-Jackson  second-sum  process.  A sixth-order  Butcher  process  is  used  to  set 
up  the  difference  table  required  before  the  second-sum  procedure  can  be  started. 


PERTURBATIONS,  TIME  AND  COORDINATE  SYSTEMS 


The  perturbing  accelerations  automatically  included  by  both  programs  are 
those  due  to  the  earth's  gravitational  potential.  The  earth's  disturbing 
function  includes  zonal  harmonics  up  to  J-  and  tesseral  harmonics  up  to  J . 
The  accelerations  due  to  atmospheric  drag  and  the  gravitational  attractions  of 
the  sun  and  moon  may  be  included  if  required.  For  a synchronous  transfer  orbit, 
the  former  is  normally  included  and  the  latter  omitted.  Atmospheric  density  is 


evaluated  using  simple  analytic  formulae  with  values  of  exospheric  temperature 
based  on  the  Jacchia  1965  model^. 


Calendar  dates  are  reckoned  in  Modified  Julian  days  (MJD) , which  are 
related  to  (ordinary)  Julian  days  by  the  formula. 


MJD  **  JD  - 2400000.5  . 

8 

The  coordinate  system  used  in  both  programs  is  that  suggested  by  Kozai  . 
The  origin  0 is  at  the  earth's  centre  of  mass  and  the  0z  axis  points  to  the 
north  pole.  Ox  lies  in  the  plane  of  the  true  equator  of  date,  but  instead  of 
pointing  to  the  true  equinox  of  date,  it  is  directed  towards  a projection  of  the 
mean  equinox  of  the  epoch  1950.0  (MJD  33281.9234);  0y  completes  the  right- 
handed  system  Oxyz  . 

5 DESCRIPTION  OF  SYNMAP 

An  overall  flowchart  for  SYNMAP  is  shown  in  Fig  3.  The  program  starts  by 
reading  all  the  input  data,  carrying  out  any  necessary  transformations  and 
setting  the  program  constants  which  are  a function  of  the  data.  The  two  random 
normal  deviate  generators  are  called  to  set  constants  from  data  stored  in  the 
program.  This  ensures  that  the  same  random  number  sequences  are  produced  in 
consecutive  runs,  thus  enabling  results  to  be  meaningfully  compared.  One 
generator  provides  a pair  of  random  normal  deviates  for  simulating  ABM  burn 
errors.  The  other  produces  six  deviates  which  are  used  to  generate  a random 
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vector  from  the  multivariate  normal  distribution  defined  by  the  tracking  error 
covariance  matrix. 

The  epoch  to  which  the  random  orbits  relate  is  stored  in  the  disc  file 
holding  the  orbits.  It  is  read  and  the  stochastic  simulation  commenced  by  read- 
ing the  first  orbit  to  be  considered.  (A  facility  exists  whereby  a simulation 
may  be  carried  out  using  every  orbit,  every  i th  orbit,  or  just  one  specified 
orbit.) 

If  tracking  errors  are  to  be  included,  a random  vector  is  generated  from 
the  error  covariance  matrix,  using  the  procedure  of  section  2.  This  vector  is 
added  to  the  orbit's  geocentric  position  and  velocity  components  (POSVEL)  at 
epoch. 


The  orbit  is  integrated  forward  until  the  ABM  burn  point  is  reached, 
ie  the  point  at  which  the  transfer  orbit  and  required  drift  orbit  plans  inter- 

f f A # # A 

sect.  At  this  point,  £ . h » 0,  where  £ is  the  radius  vector  and  h a unit 
vector  normal  to  the  required  drift  orbit  plane.  The  integration  steps  on  either 
side  of  this  point  are  identified  and  an  interpolation  performed  to  find  the 
burn  time  and  the  transfer  orbit  POSVEL.  A subroutine  is  called  to  find  the 
nominal  ABM  firing  direction  for  the  strategy  to  be  employed.  The  alternatives 
available  are  discussed  in  section  7. 

The  right-ascension  (6)  and  declination  (<J>)  of  the  nominal  firing  direc- 
tion, denoted  by  the  unit  vector  £ in  Fig  4,  are  found.  The  actual  firing 
direction  is  obtained  by  adding  errors  selected  from  a circular  distribution. 
These  errors,  e and  3 , are  shown  in  Fig  4 and  defined  as  follows.  The  actual 
firing  direction  is  assumed  to  lie  on  the  surface  of  a cone  of  half-angle  c 
whose  axis  is  the  nominal  firing  direction.  The  angle  6 is  the  azimuth  angle, 
measured  from  true  north,  and  lies  in  the  range  0 to  2ir.  Since  all  values  of  3 
are  equally  likely,  each  one  is  found  by  generating  a random  number  in  the  range 
0 to  1 and  scaling  by  2ir. 
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2 2 

2o  ■ - Eq  In  (n/m) 

Individual  values  of  e are  found  from  the  equation 

e * (j-  2o2  In  (1  - r)]^ 

where  r is  a random  variable  selected  from  a uniform  distribution  defined  on  the 
the  interval  (0,1). 

The  actual  firing  direction  is  assumed  to  have  direction  cosine 
(x,y,z)  * (1,0,0)  in  the  axis  system  whose  x-axis  lies  along  the  actual  firing 
direction,  whose  y-axis  is  in  the  plane  containing  the  nominal  and  actual  firing 
directions  and  whose  z-axis  completes  a right  handed  set.  The  direction  cosines 
^x4’y4,Z4^  actua*  firing  direction  in  the  inertial  axis  system,  are 

obtained  by  performing  a series  of  four  rotations  given  by  the  equation: 


X1 

X 

cos 

(~e)  + y sin  (-e) 

9 

yl 

- 

y 

cos 

(~e)  - x sin  (-e) 

* 

Z1 

as 

z 

• 

X2 

s 

X1 

, 

y2 

8 

yl 

cos 

(B  “ tt/2)  + Zj  sin 

(B  ~ tt/2)  , 

z2 

8 

Z1 

cos 

(B  “ tt/2)  - y(  sin 

(B  “ -rr/2)  , 

X3 

m 

X2 

cos 

(tt/2  - $)  - y sin 

(x/2  - *)  , 

y3 

8 

y2 

• 

Z3 

8 

Z2 

cos 

(ir/2  - $)  + x2  sin 

(tt/2  - <J>)  , 

■ x3  cos  (-9)  + y^  sin  (-0)  , 


y^  - y3  cos  (-9)  - x3  sin  (-0)  , 

and 
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The  original  unperturbed  orbit  is  now  integrated  forward  to  the  burn  time. 
The  actual  ABM  velocity  increment,  6v'  , is  given  by 

6v'  = 6v  + zo. 

6v 

where  z is  a random  normal  deviate  and  o.  is  the  standard  deviation  of  the 

<Sv 

nominal  ABM  velocity  increment.  The  satellite's  drift  orbit  components  after 
the  burn  are  given  by 

v,  = v.  + 6v'e 
-d  — b — 

where  v^  and  v^  are  the  drift  and  transfer  orbit  velocity  vectors  and  £ 
is  a unit  vector  in  the  actual  ABM  firing  direction. 

A subroutine  is  called  to  find  the  delta-velocity  required  for  station 
acquisition.  The  strategy  to  be  employed  will  be  mission  dependent.  Two  poss- 
ible strategies,  available  in  the  program,  are  described  in  section  8. 

Finally,  information  required  for  statistical  output  is  stored  and  the  next 
simulation  begun. 

6 TRACKING  ERRORS 

Tracking  errors  may  be  included  in  the  simulation  provided  that  an  error 
covariance  matrix  is  available  at  epoch.  This  could  be  obtained  from  a tracking 
analysis  program  using  observations  containing  errors  relating  to  the  actual 
mission. 

The  covariance  matrix  is  decomposed  into  the  product  of  a lower  triangular 
• . T 

matrix,  C , and  its  transpose,  C . For  each  simulation  the  procedure  of 
section  2 is  followed.  An  error  vector  Gz_  is  generated  and  added  to  the 
position  and  velocity  components  to  give  the  actual  orbit  at  epoch. 

7 APOGEE  MOTOR  FIRING  STRATEGY 

The  program  provides  a choice  of  three  apogee  motor  firing  strategies. 

The  ABM  may  be  used  to  give  either  a specified  drift  rate  or  a specified  flight 
path  angle  immediately  after  burn.  Alternatively  a procedure  may  be  used  to 
optimise  the  flight  path  angle  so  that  station  acquisition  can  be  achieved  with 
a minimum  total  velocity  increment.  The  required  drift  orbit  is  specified  in 
terms  of  its  inclination  (id)  and  right-ascension  («d>- 


7.1  Fixed  drift  rate 


The  drift  orbit  semi  major  axis  (a),  corresponding  to  the  required  drift 
rate  (n) , is  given  by  the  equation 

4 


aQ/(l  + n/nQ)3 


where  a^  is  the  synchronous  orbit  radius  and  n^  is  the  earth's  rotation  rate. 


The  drift  orbit  velocity  (v.)  is  given  by 

a 


vd  " 0<2/rb " 1/a)3* 


where  p is  the  earth's  gravitational  constant  and  r is  the  satellite's 

D 


radial  distance  at  ABM  burn.  Let  the  direction  cosines  of  v.  and  h , a unit 

— a — 


vector  normal  to  the  drift  orbit  plane,  be  (£  , i , 1 ) and  (h  , h , h ) 

x y z x y z 

respectively.  Then  (h  , h , h ) ■ (sin  i,  sin  fi  , - sin  i cos  ft, , cos  i,). 

xyz  ad  dd  d 

If  the  direction  cosines  of  the  transfer  orbit  velocity  vector  v^  at  the  ABM 


burn  point  are  (b^,  b^,  b^),  then 


h *.„  + h *■  + h Jl 

X X y y z Z 


+ l2  + l2 
xyz 


0 , 

1 . 


2 2 2 

and  cos  0 «■  b Jl  ♦ b l + b 1 ■ (v,  + v,  - 6v  )/2v.v,  where  0 is  the  angle 

xxyyzz  bd  bd 


between  v,  and  v,  . 
-b  -d 


Therefore 

and 


l --fhfc+hJl  1/h  , 

z l x x y y z 


[hi  + h l ] 

b(  + b l - b — 2-S *-£- 

xx  y y z h 


cos  0 


ve 


[*.  - ¥4.  * [\  - ¥4 


cos  0 


or 


Afc  + Bfc  ■ cos  0 . 

x y 


Setting  ■ [cos  0 - Bfl,y]/A  , we  have 

*1  + *?.  + fhV  + 2hiU  + *2h2l/h?  - 1 . 
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- 1 *j*[*  * s * - * *j*,  * *i[>  * a 


CSL  + DJ,  + E - 0 . 

y y 

This  equation  has  two  possible  solutions.  These  are  evaluated  and  the 
corresponding  values  of  i and  l computed.  The  drift  orbit  eccentricity 

x y A 

for  each  solution  is  calculated  and  the  unit  vector  £ in  the  ABM  firing 
direction  is  determined  for  the  solution  giving  the  minimum  drift  orbit  eccen- 
tricity using  the  equation, 

A 

<5v  s = v,  - v,  , 

where  6v  is  the  nominal  ABM  velocity  increment. 

7.2  Fixed  flight  path  angle 

Fig  5 shows  the  vectors  and  angles  referred  to  below.  The  unit  vector  u 
is  defined  by  the  equation 

u - r x h 
”b  — 

where  rfa  is  a unit  vector  along  the  radius  vector  at  the  ABM  burn  point.  The 

flight  path  angle,  y , is  measured  from  u in  the  plane  defined  by  u and 

If  v,  is  a unit  vector  along  the  drift  orbit  velocity  vector  v.  , then 
-d  “d 


u cos  y ♦ r,  sin  y 


|vdl  - |vb  • yd|  ♦ [fiv2  - (V2  - |vb  . vd|2)J*  , 


£ being  found  as  above. 


7.3  Optimised  flight  path  angle 


For  this  strategy,  the  drift  orbit  flight  path  angle  is  optimised  so  that 
the  total  velocity  increment  required  for  station  acquisition  is  a minimum. 

The  angle  Aa  , between  the  transfer  and  drift  orbit  planes  is  given  by 


> 
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cos  A ■ sin  i.  sin  i cos  An  + cos  i cos  i , 
a d t d t 


where  An  is  Che  required  nodal  rotation  and  i^  is  the  transfer  orbit  inclina- 
tion. The  angle  between  the  transfer  orbit  velocity  and  radius  vectors  is 

given  by 

cos  *t  - Yb  . £b  . 

The  angle  0 between  the  drift  orbit  and  transfer  orbit  velocity  vectors 
is  given  by 


C08  0 * COS  C08 


(>d  + sin  sin  $d  cos  Aa 


where  is  the  angle  between  the  drift  orbit  radius  and  velocity  vectors. 

2 2 2 

Now,  6v  - vd  + “ 2vdVb  C°S  8 and  80  ’ 8*ven  0 * vd  can  be  determi-ned* 

The  direction  cosines  of  v.  , (&  , l , £ ),  are  obtained  using  the  procedure 

— u x y z 

described  in  section  7.1. 


For  each  alternative  solution,  the  product  v 

“d 


v is  formed  and  the 
— b 


components  of  v^  found  for  the  scalar  product  which  is  closest  in  value  to 
cos  9 . Sufficient  information  is  available  to  enable  the  total  delta-velocity 
required  for  station  acquisition  to  be  calculated. 

The  program  sets  an  initial  value  of  <j>,  , (80°),  and  determines  the 

d 

velocity  increment  required  for  station  acquisition.  The  angle  <f>d  is  then 

incremented  in  small  steps  until  the  velocity  increment  passes  through  e minimum. 

The  step  size  is  then  reduced  and  the  procedure  repeated  until  $ is  known  to 

d 

an  acceptable  accuracy.  The  appropriate  set  of  direction  cosines  are  found  as 
described  above,  and  the  nominal  ABM  firing  direction  derived  from  them. 

8 STATION  ACQUISITION  STRATEGY 

The  station  acquisition  strategy  is  dependent  on  the  mission  requirements. 
For  most  missions,  a special  subroutine  will  have  to  be  written  and  incorporated 
in  the  program.  However,  two  alternative  strategies  are  included  and  either  may 
form  the  basis  of  a new  subroutine.  These  strategies  were  the  ones  required  for 
the  SKYNET  3 1 and  MAROTS2  missions. 

The  main  aim  of  station  acquisition  is  to  achieve  a circular  synchronous 
orbit  with  the  satellite  at  a given  longitude.  This  must  be  done  within  a speci- 
fied time  interval  and  with  the  minimum  expenditure  of  fuel  consistent  with  any 
other  constraints.  For  SKYNET  3 there  were  no  other  constraints  and  so  the  total 
velocity  increment  required  for  station  acquisition  was  minimised.  For  MAROTS 


A 


however,  it  was  required  that  the  satellite  should  acquire  station  by  the 
shortest  path.  Possible  constraints  for  other  missions  are:  (a)  moving  a 
satellite  so  that  it  is  always  visible  from  a given  ground  station;  (b)  acquiring 
station  in  one  direction  only  and  (c)  varying  the  drift  rate  during  the  station 
acquisition  sequence. 

Whichever  strategy  is  required,  the  first  steps  are  to  find  the  osculating 
elements  of  the  drift  orbit,  the  satellite's  longitude  at  ABM  burn,  its  longitude 
after  any  tracking  period  has  elapsed  and  its  drift  rate  relative  to  the  earth. 

In  general,  the  ABM  firing  strategy  will  have  been  chosen  to  minimise  the  drift 
orbit  eccentricity,  although  a correction  may  still  be  required.  A manoeuvre 
may  also  be  required  to  rotate  the  orbit  plane.  This  is  especially  likely  for 
missions  which  require  north-south  stationkeeping. 

All  delta-velocity  calculations  are  made  with  linearised  equations. 

8. 1 SKYNET  3 station  acquisition 

The  velocity  increments  needed  to  acquire  station  using  both  easterly  and 
westerly  drifts  are  determined  and  the  direction  giving  the  minimum  delta- 
velocity  selected.  If  the  initial  drift  rate  is  inadequate  for  the  satellite  to 
reach  the  specified  longitude  in  the  prescribed  time,  a manoeuvre  is  performed. 
The  velocity  increment  required  to  stop  the  satellite  ic  then  computed.  It  is 
assumed  that  these  manoeuvres  are  also  used  to  reduce  the  orbital  eccentricity. 
Finally,  the  velocity  increments  required  to  remove  any  residual  eccentricity 
and  perform  an  inclination  correction  are  calculated. 

8.2  MAROTS  station  acquisition 

The  drift  direction  which  gives  the  shortest  path  to  station  is  identified 
and  the  satellite's  initial  drift  rate  examined.  If  it  is  inadequate,  or  in  the 
wrong  direction,  a manoeuvre  is  performed.  A further  manoeuvre  is  made  to  stop 
the  satellite  on  station.  Again  it  is  assumed  that  these  manoeuvres  are  also 
used  to  reduce  the  orbital  eccentricity.  Lastly,  the  velocity  increments 
required  to  remove  any  residual  eccentricity  and  perform  an  inclination  correc- 
tion are  determined. 

9 STATISTICS 

The  program  computes  the  means  and  standard  deviations  of  the  right- 
ascension  and  declination  of  the  required  ABM  firing  direction;  the  solar  aspect 
angle;  the  spin  axis  elevation  angle  (defined  as  the  angle  between  the  actual  ABM 
firing  direction  and  the  plane  normal  to  the  radius  vector  at  the  ABM  burn  point) 
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the  drift  orbit  osculating  elements;  the  longitude  of  the  burn  point  and  the 
total  velocity  increment  (AV)  required  for  station  acquisition.  The  following 
expressions  are  used  to  evaluate  means  and  standard  deviations 


i-n 


where  n is  the  number  of  simulations.  A count  is  also  kept  of  the  number  of 
violations  of  the  solar  aspect  angle  constraint. 

For  representative  results  to  be  obtained,  at  least  250  simulations  should 
be  performed.  For  more  detailed  analysis,  it  is  suggested  that  at  least  1000 
cases  be  considered. 

Since  the  distributions  of  some  quantities,  notably  AV  , are  skew,  no 
great  significance  should  be  attached  to  the  actual  values  of  mean  plus  three 
standard  deviations,  ie  x + 3a^  . (For  a normal  distribution,  the  probability 
of  a quantity  not  exceeding  this  value  is  0.997.) 

The  program  also  produces  a histogram  of  the  AV s required. 

10  INPUT  DATA  FOR  P0INT2 

10.1  Overall  description 

For  convenience,  the  input  data  is  treated  here  as  a set  of  17  punched 
cards.  The  order  of  the  cards  is  as  follows: 

(i)  time  card 

(ii)  control  cards 

(iii)  nominal  orbit  card 

(iv)  gross  attitude  manoeuvre  card 

(v)  covariance  matrix  or  equivalent  lower  triangular  matrix. 


All  records,  except  those  specifying  a matrix,  are  read  in  free  format. 
In  the  description  which  follows,  parameters  and  constants  are  referred  to  by 
their  Fortran-variable  names.  A typical  data  deck  is  illustrated  in  Fig  6. 


i 


10.2  Time  card 


The  time  card  contains  four  parameters  specifying  the  injection  epoch  and 

the  times  at  which  the  sets  of  random  orbits  are  required. 

MJDOCH  The  MJD  number  of  the  epoch  at  which  injection  into  the  transfer  orbit 
takes  place. 

EP  The  time  of  injection,  expressed  as  a fraction  of  a day  relative  to 

MJDOCH. 

TINT1  The  time  at  which  the  first  set  of  random  orbits  is  required  in  hours 
from  epoch. 

TINT2  The  time  at  which  the  second  set  of  random  orbits  is  required,  in  hours 
from  epoch.  TINT2  ■ 0.0  if  only  one  set  of  orbits  is  required. 

10.3  First  control  card 

This  card  contains  three  parameters  (all  integer)  which  control  the  work- 
ing of  the  program. 

NC  The  number  of  orbits  required  in  each  set. 

NOLOT  If  NOLOT  = 2,  the  program  reads  an  error  covariance  matrix  correspond- 

ing to  a set  of  injection  conditions,  with  units  ft,  ft/s  and  degrees. 
The  units  are  converted  to  km,  km/s  and  radians  and  the  matrix  decom- 
posed into  its  lower  triangular  form.  If  NOLOT  ■ 1 , the  lower  triangu- 
lar matrix  is  read  directly. 

LSP  If  LSP  * 1,  luni-soiar  perturbations  are  included  in  the  integration. 

10.4  Second  control  card 

This  card  contains  one  control  parameter,  one  parameter  which  is  used  as 

both  a data  item  and  a control  parameter,  and  two  data  parameters. 

HN  Although  the  integration  step-length  (H)  is  set  in  the  BLOCK  DATA 

segment,  it  may  be  modified  to  H/HN  by  setting  HN  0. 

2 

ARMA  The  area-to-mass  ratio  of  the  satellite  (m  /kg).  If  set  to  0.0,  air- 
drag  terms  are  excluded  from  the  integration. 

SMEAN  The  value  of  the  solar  10.7cm  radiation  flux  averaged  over  three  days 

-22  -2  -1 

and  measured  in  units  of  10  Vta  Hz  . If  ARMA  - 0.0,  then 
SMEAN  - 0.0. 
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SOLBAR  The  value  of  the  solar  10.7  cm  radiation  flux  averaged  over  three  months 

-22  ~2  - l 

and  measured  in  units  of  10  Wm  Hz  If  ARMA  ■ 0.0,  then 
SOLBAR  =»  0.0. 

10.5  Nominal  orbit  card 


i 


j 

* 


The  nominal  orbit  at  epoch  is  supplied  as  a set  of  injection  conditions: 

IPN(6)  Velocity;  climb  angle;  azimuth;  radial  distance;  latitude  and  longitude 
of  injection,  in  units  km,  km/s  and  degrees. 

10.6  Gross  attitude  manoeuvre  card 

This  card  contains  one  parameter  which  is  used  both  as  a control  parameter 
and  a data  item,  and  three  data  parameters. 

NAM  The  apogee  number  at  which  the  gross  attitude  manoeuvre  is  performed. 

If  NAM  ■ 0,  no  manoeuvre  is  included. 

DVT  The  transverse  delta-velocity  component  resulting  from  the  manoeuvre 

(km/s) . 

DVN  The  normal  delta-velocity  component  resulting  from  the  manoeuvre  (km/s). 

DVR  The  radial  delta-velocity  component  resulting  from  the  manoeuvre  (km/s). 

10.7  Covariance  or  lower  triangular  matrix 

The  matrix  is  supplied  on  12  cards,  each  containing  three  elements  punched 
in  3E15.8  format.  Cards  one  and  two  hold  the  first  row,  cards  three  and  four 
the  second  row  and  so  on. 

1 1 INPUT  DATA  FOR  SYNMAP 

The  input  data  is  divided  into  three  categories.  The  first  is  a card  data 
deck  containing  information  about  the  mission.  The  second  is  a set  of  random 
orbits  stored  in  an  unformatted  disc  file.  The  third  is  data  set  within  the 
program.  A typical  card  deck  is  illustrated  in  Fig  7. 

1 1 . 1 Data  deck 

The  data  deck  contains  a minimum  of  four  cards,  plus  an  additional  12 
cards  if  tracking  errors  are  to  be  included.  The  order  of  the  cards  is  as 
follows:- 


l 
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(i)  control  card 

(ii)  apogee  manoeuvre  card 

(iii)  station  acquisition  card 

(iv)  perturbation  card 

(v)  tracking  covariance  matrix. 

All  records  are  punched  in  free  format  except  where  indicated.  Parameters 
and  constants  are  referred  to  by  their  Fortran  variable  names. 

11.1.1  Control  card 

This  card  contains  five  parameters  which  control  the  working  of  the 
program: 

NC  Number  of  cases  (ie  simulations). 

STR  If  STR  ■ 0.0,  ABM  pointing  errors  are  neglected.  Otherwise  STR  is  the 

starting  value  for  the  RAND0MN0  function  and  is  a positive  odd  integer. 

NORB  If  NORB  1,  all  random  orbits  are  used.  If  NORB  > 1 and  NC  = 1 , only 
orbit  number  NORB  is  used.  If  NC  * i and  NORB  * j,  then  orbits 
j,  2 j , 3 j , ....  ij  are  used. 

NM  If  NM  0,  a tracking  data  covariance  matrix  is  read  and  tracking 

errors  are  included. 

NSTATS  If  NSTATS  > 0,  statistical  data  is  output. 

11.1.2  Apogee  manoeuvre  card 

This  card  contains  the  following  variables.  The  first  is  punched  in  A4 
format  and  the  others  in  free  format. 

SUB 

DELTAV 
DIN 
MIA 
GAM 


The  name  of  the  subroutine  to  be  used  to  determine  the  ABM  firing 
direction. 

The  nominal  ABM  velocity  increment  (km  s *). 

The  required  drift  orbit  inclination  (degree). 

The  required  drift  orbit  right  ascension  (degree)  . 

If  SUB  ■ DR,  GAM  is  the  required  drift  orbit  drift  rate  (degree/s). 

If  SUB  s FPA,  GAM  is  the  required  drift  orbit  flight  path  angle  (degree). 
If  SUB  ■ OFPA,  GAM  is  not  used. 
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11.1.3  Station  acquisition  card 

This  card  contains  the  following  variables.  The  first  is  punched  in  A4 
format  and  the  others  in  free  format. 

SUBS  The  name  of  the  station  acquisition  subroutine  to  be  used. 

TAC  The  total  time  available  for  station  acquisition  (days). 

TR  Tracking  time  required  after  ABM  burn  (days). 

PHI  Required  station  longitude  (degrees  east  of  Greenwich). 

XIL  I Minimum  and  maximum  permitted  values  of  drift  orbit  inclination 
XIH  j (degrees) . 

11.1.4  Perturbation  card 

This  card  contains  the  following  variables: 

Cl  Although  the  integration  step-length,  H,  is  set  in  the  BLOCK  DATA  seg- 
ment, it  may  be  modified  to  H = H/Cl  by  setting  Cl  # 0. 

LSP  If  LSP  > 0,  luni-solar  perturbations  are  included  in  the  integration. 

2 

ARMA  The  area-to-mass  ratio  of  the  spacecraft  (m  /kg).  If  ARMA  > 0,  air 
drag  terms  are  included  in  the  integration. 

SMEAN  The  value  of  the  solar  10.7  cm  radiation  flux  averaged  over  three  days 

-22  -2  -1 

and  measured  in  units  of  10  W m Hz  If  ARMA  ■ 0.0,  then 
SMEAN  - 0.0. 

SOLBAR  The  value  of  the  solar  10.7  cm  radiation  flux  averaged  over  three 

. . -22  -2  -I 

months  and  measured  in  units  of  10  W m Hz  . If  ARMA  =*  0.0,  then 
SOLBAR  - 0.0. 

11.1.5  Tracking  covariance  matrix 

If  tracking  errors  are  to  be  included,  this  matrix  is  supplied  on  12  cards 
punched  in  3E0.0  format.  Cards  one  and  two  hold  row  one,  cards  three  and  four 
row  two  and  so  on. 

11.2  Random  orbits 

The  epoch  and  random  orbits  are  supplied  in  an  unformatted  disc  file.  The 
epoch  is  specified  as  a Modified  Julian  day  number  (MJD)  and  a fraction  of  a day 
relative  to  it.  Each  orbit  is  specified  as  a set  of  geocentric  position  and 
velocity  components  (x,  y,  z,  x,  y,  z,)  in  km  and  km  s * . 
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11.3  Data  stored  in  the  program 

Four  data  items  are  set  in  the  program  itself.  The  BLOCK  DATA  segment  sets 

the  maximum  and  minimum  permitted  values  of  solar  aspect  angle  (115°  and  65° 

respectively).  Data  statements  in  the  MASTER  set  values  of  (in  degrees)  and 

the  standard  deviation  of  the  ABM  velocity  increment,  is  defined  as  the 

value  of  £ (see  section  5)  which  is  only  exceeded  in  three  cases  out  of  1000. 

2 

The  value  of  2o  is  deduced  from  by  setting 

0.997  ■ 1 - exp(-  e^/2o2)  , 

ie 

2 

2a2  - 

a In  (0.003) 

The  standard  deviation  of  the  ABM  velocity  increment  is  divided  by  the  nominal 
velocity  increment  to  give  a value  in  the  range  0 to  1 . 

12  SUN-MOON  DATA  FOR  BOTH  PROGRAMS 

To  enable  the  solar  aspect  angle  to  be  calculated  and  luni-solar  perturba- 
tions to  be  included  in  the  integration,  a table  of  sun-moon  coordinates,  at 
daily  intervals  with  second  and  fourth  differences,  must  be  supplied  in  a disc 
file.  A fuller  description  is  given  in  the  specification  of  subroutine  SMPOS 
in  Appendix  B. 

13  DESCRIPTION  OF  SYNMAP  OUTPUT 

A sample  line  printer  output  is  shown  in  Fig  8.  A heading  is  first  printed, 
giving  information  about  the  case  being  considered.  A minimum  of  two  lines  is 
then  output  for  each  simulation.  If  the  solar  aspect  angle  constraint  has  been 
violated,  a warning  is  printed  to  this  effect.  The  uext  line  consists  of  the 
case  number,  right  ascension  and  declination  of  the  actual  ABM  firing  direction, 
the  drift  orbit  osculating  elements,  the  longitude  of  the  burn  point  and  the 
satellite's  drift  rate.  This  is  followed  by  one  line  for  each  drift  direction 
considered  for  station  acquisition.  The  direction  is  printed  together  with  the 
selected  drift  rate,  the  delta-velocities  required  to  set  up  the  drift  rate  and 
stop  the  satellite  on  station,  the  sum  of  the  last  two  quantities,  the  delta- 
velocity  needed  to  circularize  the  orbit,  the  delta-velocity  required  for 
inclination  correction  and  the  total  delta-velocity  required  for  station 
acquisition.  Finally,  when  the  simulation  is  complete,  the  statistical  informa- 
tion and  the  histogram  of  delta-velocities  is  output. 
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Appendix  A 
POINT2  PROGRAM  UNITS 


reduces  an  angle  to  the  range  0 to  2ir 
reduces  an  angle  to  the  range  -it  to  it 
generates  random  ndftnal  deviates 
sets  certain  constants 

integrates  second  order  differential  equations 
simulates  gross  attitude  manoeuvre 
finds  a named  area  on  a disc  file 

converts  a number  into  its  integer  and  fractional  parts 

converts  injection  parameters  to  coordinates 

current  time  in  seconds  from  midnight 

performs  triangular  decomposition  of  a matrix 

computes  atmospheric  density 

controls  orbit  integration  and  interpolation 

auxiliary  for  DEQRSPT,  computes  perturbing  accelerations 

forms  scalar  product 

reads  sun/jpoon  coordinates  from  disc  and  interpolates 
performs  interpolation 

determinA  polar  coordinates  from  cartesians  (two-dimensional) 
permits  disc/core  transfers. 


The  subprograms  INFIND,  INTFRC  and  SCPROD  are  written  in  PLAN  and  are  used 
as  semi-compiled  segments. 


(b) 


PDAUXP  appears  in  the  subprogram  (DEQRSPT)  which  calls  it  as  DAUX.  This 
is  organist  through  subprogram  arguments  with  an  EXTERNAL  statement  in 
subroutine  0RBIT2.  The  object  is  to  permit  substitution  of  a different 
subroutine  without  alteration  of  the  calling  program. 


(c) 


The  subprograms  ITIME  and  UTD4,  though  not  part  of  standard  FORTRAN,  are 
provided^ptomatically  by  the  1900  series  FORTRAN  compilers  and  are  not 


described  here. 
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Summary 

Language 

Author 

Function  statement 
Input  argument 
X 

Output  function 
ANGLE 

Use  of  COMMON 

Source  deck 

Local  storage  used 

Subordinate  subprograms 

Explanation 

fractional  part  of  x/2ir 


FUNCTION  ANGLE 

The  function  reduces  an  angle  x (in  radians)  to  the 
range  0 < x < it  . 

1900  Fortran. 

Diana  W.  Scott  (April  1969). 

FUNCTION  ANGLE (X). 

Angle  x in  radians. 

Value  of  x ± 2mr  , such  that  0 ANGLE  < 2tt  . 

None. 

6 cards,  including  1 comment  card  (ICL  code). 

1 real  variable. 

None. 

The  standard  function  AMOD  is  used  to  give  the 
If  this  is  negative,  2n  is  added  to  the  result. 
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Appendix  B 
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FUNCTION  ANGLI 


Summary 

- The  function  reduces  an  angle  x to  the  range 

-it  < x < n . 

Language 

- 1900  Fortran. 

Author 

- Diana  W.  Scott  (April 

1969). 

Function  statement 

- FUNCTION  ANGLI (X) 

Input  argument 

- 

X 

Angle  x in  radians. 

Output  function 

- 

ANGLI 

The  value  of  x ± 2nir 

such  that  -it  < ANGLI  < it 

Use  of  COMMON 

- None. 

Source  deck 

- 7 cards,  including  one 

common  card  (ICL  code). 

Local  storage  used 

- 1 real  variable. 

Subordinate  subprograms 

- None. 

Explanation  - The  standard  function  AMOD  is  used  to  give  the 

fractional  part  of  x/ 2tt  . If  this  is  greater  than  n , 2ir  is  subtracted, 
it  is  less  than,  or  equal  to,  -it  , 2ir  is  added. 


If 
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SUBROUTINE  ARANOR 


Summary 


Language 


Author 


Subroutine  statement 


~ The  subroutine  generates  pseudo-random  numbers 
normally  distributed  with  mean  zero  and  standard 
deviation  one. 

- 1900  Fortran 

- A.W.  Odell  (November  1973), 

- SUBROUTINE  ARANOR  (Z,  N) . 


Input  argument 


Use  of  COMMON 


Integer  specifying  mode  of  operation  and/or  number 
of  random  numbers;  see  explanation. 


Input  and  output  argument  - 


Array  containing  random  numbers;  see  explanation. 


- None. 


Source  deck 


- 19  cards  (ICL  code). 


Subordinate  subprograms 


Explanation 


Local  storage  used  - 3 real  variables  and  1 integer  variable. 

Subordinate  subprograms  - The  subroutine  ITIME  (standard  ICL  1900  library 

subroutine) . 

Explanation  - The  method  used  is  described  in  Ref  5.  Two 

sequences  of  rectangularly  distributed  random  numbers  are  generated  using  the 
equations 

x.  , = 47 lOlx.  + 1410065412  (mod  231  - 1) 

l+l  i ' 

yi+l  * 46893Fi  + 2115099118  (mod  231  - 1). 

A sequence  of  normally  distributed  random  numbers  is  obtained  from 

z.  - (-2  log  (x./(231  - 1))J  sin  (2ny./(231  - l)). 

To  obtain  a set  of  n random  numbers,  N may  be  set  to  n before  entry.  On 
exit  from  the  routine,  Z will  contain  the  random  numbers.  The  initial  x and  y 
will  have  been  obtained  using  the  computer  clock  and  the  set  of  numbers  will  not 
therefore  be  repeatable. 


i 
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In  order  to  provide  a repeatable  set  of  numbers,  or  if  the  subroutine  is  to 
be  run  on  a computer  without  a clock,  x and  y should  be  initialised  by  placing 
suitable  values  in  Z(l)  and  Z(2)  and  calling  the  subroutine  with  N - 0.  A set  of 
n random  numbers  can  now  be  obtained  in  Z by  re-entering  the  subroutine  with  N set 
to  -n.  The  use  of  the  subroutine  with  N negative  may  also  be  used  to  continue  a 
set  of  numbers  produced  by  the  last  entry  to  the  subroutine. 
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BLOCK  DATA  (POINT) 

Simmary  - The  Fortran  BLOCK  DATA  segment  is  used  to  set  initial 

values  for  certain  quantities  stored  in  conmon  blocks 
/PETURB/,  /CINTEG/,  /CONST/  and  /CON/. 


Language  - 1900  Fortran. 


Author 

- M.D.  Palmer  (January  1975). 

Data  in  /PETURB/ 

- 

Variable  name 

Value 

Explanation 

HALF CD 

1.1 

Product  JCD  . 

RATE 

1.1 

Ratio  of  angular  velocity  of  the  atmosphere 

to  that  of  the  earth. 

Data  in  /CINTEG/ 

- 

Variable  name 

Value 

Explanation 

H 

2s/96 

Integration  step  length. 

PD  (18) 

All  elements  0.0 

Partial  derivatives  of  position. 

PDVEL  (18) 

All  elements  0.0 

Partial  derivatives  of  velocity. 

Data  in  /CONST/ 

- 

Variable  name 

Value 

Explanation 

EMU 

398601. 3kmV2 

Earth's  gravitational  constant. 

EJ2 

1082.637E-6 

Earth's  second  zonal  harmonic  J2  . 

EJ3 

-2.53I0E-6 

Earth's  third  zonal  harmonic  J^  • 

EJ4 

-I.6190E-6 

Earth's  fourth  zonal  harmonic  J.  . 

4 

C22 

2.4369E-6 

Tesseral  harmonic  coefficient  C22  • 

S22 

-I.4005E-6 

Tesseral  harmonic  coefficient  S22  . 

C33 

0. 7387E-6 

Tesseral  harmonic  coefficient  C^  • 

S33 

1.4343E-6 

Tesseral  harmonic  coefficient  . 

C44 

-0. 1846E-6 

Tesseral  harmonic  coefficient  C, . . 

44 

S44 

0.2508E-6 

Tesseral  harmonic  coefficient  S^  . 

C31 

2.0I92E-6 

Tesseral  harmonic  coefficient  . 

S31 

0. 2278E-6 

Tesseral  harmonic  coefficient  S^  . 

C42 

0. 3444E-6 

Tesseral  harmonic  coefficient  • 

S42 

0.7021E-6 

Tesseral  harmonic  coefficient  • 

ERAD 

6378.163km 

Earth's  mean  equatorial  radius. 
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Variable  name 

Value 

EOMEGA 

7.2921 151 47E-5rad/ 

PREC 

6.079E-12rad/s 

SUNMU 

1 . 3271 27E1 lkm3s~2 

SELMU 

4902. 756km3s~2 

Data  in  /CON/ 

- 

Variable  name 

Value 

EJ5 

-0.246E-6 

EJ6 

0. 558E-6 

EJ7 

-0.326E-6 

EJ8 

-0.209E-6 

EJ9 

-0.094E-6 

C43 

1 .0390E-6 

S43 

-0. 1 192E-6 

C41 

-0.5I75E-6 

S41 

-0.481 4E-6 

C32 

0. 7783E-6 

S32 

-0. 7552E-6 

Source  deck  - 18 

cards  (ICL  code). 

Explanation 

Earth's  angular  velocity. 

Earth's  precession  rate. 

Sun's  gravitational  constant. 
Moon's  gravitational  constant. 

Explanation 

Earth's  fifth  zonal  harmonic  Jj  . 

Earth's  sixth  zonal  harmonic  Jg  . 

Earth's  seventh  zonal  harmonic 

Earth's  eighth  zonal  harmonic  Jg . 

Earth '8  ninth  zonal  harmonic  . 

Tesseral  harmonic  coefficient  C., 

43 

Tesseral  harmonic  coefficient  S,_ 

43 

Tesseral  harmonic  coefficient  C. , 

41 

Tesseral  harmonic  coefficient  S.. 

4 1 

Tesseral  harmonic  coefficient 
Tesseral  harmonic  coefficient  S.^ 
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SUBROUTINE  DEQRSPT 

Summary  - The  subroutine  integrates  a set  of  up  to  22  simultaneous 

second-order  differential  equations  of  the  form 

y ^ m j • • • • y 1 1 • • • *y^)  > ^ * i>«»«>n  , us  mg 

a Gaussian  eighth-order  second-sum  predictor-corrector^. 
The  integration  is  started  using  a Butcher  sixth-order 
seven  stage  Runge-Kutta  process10  to  set  up  the  differ- 
ence table  required  before  the  second-sum  procedure  can 
take  over. 


Language 


- 1900  Fortran. 


Authors  - A.W.  Odell  and  G.J.  Davison  (April  1973). 

Subroutine  statement  - SUBROUTINE  DEQRSPT  (NTRY,  DAUX) . 


Input  arguments 
NTRY 

DAUX 

Output  arguments 
Use  of  COMMON 


1 for  an  initialization  entry  to  the  subroutine,  and 

2 for  all  normal  entries  (see  Explanation). 

The  auxiliary  subroutine  which  evaluates  y.  . 

- None. 

- Certain  quantities  in  common  block  /CINTEG/  are  used 
as  follows :- 


Input  arguments  in  /CINTEG/  - 

NEQ  Number  of  equations  (normally  4 or  22). 

H Integration  step  size,  h (positive  or  negative). 

Input  and  output  arguments  in  /CINTEG/  - 

T Independent  variable,  t . 

Y(22)  Dependent  variables,  y^  . 

YP(22)  First  derivatives,  y,,  . 

Y2P(22)  Second  derivatives,  y^  , as  computed  by  the  auxiliary 

subroutine  DAUX. 

Output  arguments  in  /CINTEG/  - 

IAUX  Set  to  -I  following  the  initialisation  entry  (NTRY>I), 

to  0 if  T has  been  changed  (during  a normal  entry), 
and  to  1 otherwise. 
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Source  deck  - 218  cards  (ICL  code). 

Local  storage  used  - 10  integer  variables,  1 logical  variable,  63  real 

variables  and  814  real-array  elements. 

Subordinate  subprograms  - The  auxiliary  subroutine,  noned  in  the  call  to 

DEQRSPT  which  takes  the  role  of  DAUX. 

Explanation  - DEQRSPT  integrates  a set  of  N(<22)  simultaneous 

second-order  differential  equations  of  the  form  y^  - f £ (t ,y ! , . . ,yN,y ( , . . .yN)  , 
i “ l,...,N  , using  a Gaussian  eighth-order  second-sum  predictor  corrector.  The 
integration  is  started  using  a Butcher  (6,7)  R-K  process  to  set  up  the  difference 
table  required  before  the  second-sum  procedure  takes  over. 

Prior  to  any  integration,  DEQRSPT  must  first  be  called  with  NTRY  - 1.  This 
is  the  initialization  entry  in  which  the  step  size  h'(-  jh)  is  set  for  the 
Butcher  integration  and  DAUX  is  called  to  evaluate  y.  at  t^  . All  subsequent 
entries  are  made  with  NTRY  ■*  2,  and  one  integration  step  (two,  whilst  still  in 
the  Butcher  mode)  is  performed  before  control  is  returned  to  the  calling  program. 

If  y.  . and  y.  . are  the  values  of  y.  and  y.  at  t ■ t,  the 

X g K 1 X 1 K 

formulae  used  in  the  next  Butcher  step  to  evaluate  y.  , , . y.  and 

„ 7i,k+l * 7i,k+l 

y.  . . , the  values  of  y. , y.  and  y.  at  t - t,  ♦ h’  , are: 

X y I XXX  K 


yi,k*l 


7 


> 


/ 

■ n.k  * E v*. 

s-l 


and 


yi,k+l  " fi(tk  + h,,yl,k+l,“*yN,k+l,yl,k+l'**,yN,k-H)  ’ i ■ 1 , . . • ,N 


where  k 


s-l 


*i.  * * Z • 

j-1 


a ■ I N , 
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Given  y.  ^ and  q * q is  obtained  from  DAUX,  and  the  firat  two 

Butcher  steps  of  size  h'  then  give  y.  , y.  and  y.  at  time 

i»i  i9l  # I 

tj  ■ tfl  ♦ h . The  process  is  repeated  until  y^  y^  g and  y^  g are 
obtained  after  a total  of  16  steps  (eight  calls  to  the  subroutine).  The 
Butcher  process  is  now  complete. 

On  the  ninth  call  (with  NTRY  - 2)  to  DEQRSPT  the  following  difference  table 
is  constructed  for  each  y^  = f^  . 
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ri,0 


i.l 


ri.1 


Vl,2 


7i.2 


fi,2 


Vi,3 


f • • 

i.3 


ri.3 


1 • / 

i.4 


t • / 

i,4 


-2 

•• 

i.3 

yi.4 

-2 

•• 

i.4 

7i!s 

yi.5 

-2 

•• 

i.5 

7l!e 

yi.6 

-2 

„ 

i.6 

7i!7 

yi»7 

-2 

•• 

i.7 

7m 

yi.8 

s 

s' 

00 
CM  • 

/ 

/ 

yi.9 

vi!9 

’i.5 


*1.3 


t ' * 

1,6 


'i.« 


Vi.7 


i.7 


”i.8 


v: 


i.« 


v: 


i.9 


7^ 

i.3 


V? 

i.5 


Vi.5 


i,5 


7? 

1,6 


7*  ft 

1,6 


7? 

1,6 


i»6 


i.7 


'1.7 


9 • *7 


i.8 


r • -» 
1.7 


76  „ 

i,8 


1 .8  - 
" v8 
" i.9 


i.7 


Vi.8 


7^ 

i »*» 


7^ 

i.8  ^ 


" 77 


i.9 


V6  o 
i.9 


^ 7 


i.9 


^ 7 


7? 

i.9 


i.9 


First  the  known  y,  are  differenced  to  give  the  part  of  cha  table  to  the  right 
#>  i „ | _2 

of  the  y , and  then  the  first  and  second  sums  7.  ^ and  7^  ^ ara  formed, 

using  the  equations 


_~2  . — 2 _ 2 _ _4  _ _6  m _8 

7i,3  " h yi,4  " Vi, 4 “ Vi, 5 " Vi, 6 " Vi, 7 “ Vi, 8 
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The  remaining  7 and  7 quantities  above  the  dotted  line  are  defined  on 
the  basis  that  the  difference  between  any  entry  and  the  entry  above  must  equal 
the  entry  on  the  right.  Actually  only  7^  g and  7,,  are  required  explicitly 
and  these  are  given  by 


A table  containing  the  coefficients  used  in  the  above  and  following  equations  is 
appended. 

The  ninth  normal  call  to  the  subroutine  continues  with  an  integration  step 
using  the  Gaussian  (predictor)  formulae 


which  require  only  the  quantities  inmediataly  above  the  dotted  line  in  the  table 


y^  g is  then  obtained  using  DAUX  and  the  row  of  differences  under  the 
diagonal  line  found.  This  row  is  then  used  to  obtain  corrected  values  of  y 
and  y using  the  equations 
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+ Vi.9 

1 2 
+ D V.  0 + D-V.  n 

1 1,9  2 i,9 

* Vi.9 

+ DCV?  . 

+ D,?* 

+ D-V?  „ 

+ D0?®  0 

5 i,9 

6 i,9 

7 i,9 

8 i,9, 

*i,9  ' h(VlU  * Vi, 9 * E.7i,9  * Vi.9  * Vi,9  * E47i,9 

* E57i,9  * Vi, 9 * ViN  * Vi.9)  • 

••  8 
y.  Q is  then  redetermined,  and  the  row  of  differences  out  to  V?  under 
1 , y i ,9 

the  dotted  line  recalculated.  This  row  is  then  used  to  obtain  the  final 
corrected  values  of  y.  Q and  y.  _ using  the  above  equations. 

Coefficients  for  A..  B..  C. . D. . E.  and  F. 


863  275 


12096  4032 


33953 


518400 


8183 


129600 


22809600 


3250433 


53222400 


9829 

407 

330157 

6048 

3628800 

172800 

159667200 

863 

275 

33953 

8183 

60480 

24192 

3628800 

1036800 

1070017 


3628800 


25713 


89600 
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SUBROUTINE  GATMAN 

Summary  - The  subroutine  integrates  a satellite  orbit  to  the 

apse  at  which  a gross  attitude  manoeuvre  is  to  be 
performed.  It  interpolates  to  find  the  geocentric 
cartesian  components  of  position  and  velocity  at 
the  apse  and  then  determines  the  satellite's  velocity 
components  immediately  after  the  manoeuvre. 

- 1900  Fortran. 

- M.D.  Palmer  (August  1975). 

- SUBROUTINE  GATMAN  (I,  TAP,  X,  Y,  Z,  XD,  YD,  ZD). 

The  Modified  Julian  day  number  of  the  gross  attitude 
manoeuvre . 

TAP  The  time  of  the  manoeuvre,  in  fractions  of  a day, 

relative  to  I. 

X,  Y,  Z The  satellite's  geocentric  position  components, 

r_  = (x,  y,  z)  , at  the  time  of  the  manoeuvre  (km). 

XD,  YD,  ZD  The  satellite's  geocentric  velocity  components, 

v = (x  , £ , z ) , immediately  after  the  manoeuvre, 

“TO  in  IQ  IQ 

(km/s) . 

Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CINTEG/, 

/CONST/,  /CSM00N/  and  /INTP/  are  used  as  follows 


Input  arguments  in  /CINTEG/  - 


MJDOCH 


Modified  Julian  day  number  of  epoch. 


PP(3)  Satellite's  geocentric  position  components, 

r_  “ (x,  y,  z) , initially  at  epoch  and  subsequently 
at  the  latest  integration  step  (km). 


PV(3)  Satellite's  geocentric  velocity  components, 

v = (x,  y,  z),  initially  at  epoch  and  subsequently 
at  the  latest  integration  step  (km/s). 


PA(3) 


Satellite's  geocentric  acceleration  components  at 

2 

the  latest  integration  step  (km/s  ) . 


i 
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Input  arguments  in  /CONST/  - 


EMU 


• • 5 —2 

Earth's  gravitational  constant  (km  s ) 


Input  arguments  in  /CSMOON/  - 


MJDT 
TIMET 

Input  arguments  in  /INTP/  - 
EP 

NAM 


Modified  Julian  day  number  of  the  current  time. 

The  current  time,  in  fractions  of  a day,  relative 
to  MJDT. 


DVT 

DVN 

DVR 


The  time  of  epoch,  in  fractions  of  a day,  relative 
to  MJDOCH. 

The  number  of  the  apse  at  which  the  manoeuvre  is  to 
be  carried  out. 

The  transverse,  normal  and  radial  components  of 
delta  velocity,  resulting  from  the  manoeuvre  (km/s) 


Output  arguments  in  /CINTEG/  - 

The  initial  values  of  dx/ds,  dy/ds,  dz/ds. 


XVEL 

YVEL 

ZVEL 


TVEL 

Source  deck 
Local  storage  used 

Subordinate  subprograms 


The  initial  value  of  dt/ds  (s~*) 

- 84  cards  (ICL  code). 

- 15  real  array  elements;  39  real  variables;  and 
2 integer  variables. 

- The  subroutines  DEQRSPT , MODAT,  PDAUXP,  SMPOS,  TRINV 
and  the  functions  ANGLE,  ANGL1 , INFIND,  INTFRC , 
SOLVIN,  SCPROD  and  UTD4. 


Explanation  - The  subroutine  starts  and  controls  the  integration 

of  a satellite  orbit  until  a specified  apse  is  reached.  An  interpolation  is 
then  performed  to  find  the  position  and  velocity  components  at  the  apse,  and  the 
satellite's  velocity  immediately  following  the  manoeuvre  is  determined. 

The  variables  TVEL(dt/ds),  XVEL(dx/ds) , YVEL(dy/ds)  and  ZVEL(dy/ds)  are 
set  in  terms  of  the  independent  variable  s.  The  integration  is  started  by  call- 
ing DEQRSPT  with  the  statement  CALL  DEQRSPT  (IND,  PDAUXP).  On  the  first  call, 
IND  - 1 and  on  the  second  and  subsequent  calls  IND  <■  2.  After  each  integration 


Appendix  B 


41 


step,  the  time  elapsed  from  epoch  is  calculated  in  days.  The  length  of  the  rad- 
ius vector,  r^  , and  its  rate  of  change,  r^  , are  determined.  If  r^  . < 0 , 

there  is  an  apse  between  the  last  two  integration  steps.  A counter  is  incremented 
at  each  apse  and  when  the  count  is  equal  to  the  variable  NAM,  an  interpolation  is 
performed  to  find  the  satellite's  position  r « (x,  y,  z)  and  velocity 
v = (x,  y,  z)  at  the  apse. 

The  equations  used  for  interpolation  are:- 


f(t) 


q3[d  + 3p  + 6p2)f(tj)  + hp(l  + 3p)f ' (t | ) + | p2f "(tj)] 

2 

+ P3[(1  + 3q  + 6q2)f(t2)-hq(l  + 3q)f'(t2)  +|  qV^,,)] 


for  position  components  and 

f'(t)  - 30h_lp2q2[f(t2)  “ f(tj)]  + q2( 1 + 5p) ( 1 - 3p)f ' (t j) 

+ P20  + 5q)(l  - 3q)f'(t2)  + | pq2 (2  - 5p)f”(t1) 


- \ p2q(2  - 5q)f'(t2) 


for  velocity  components,  where  ll‘'C2~ti  » P“  Ct~  tj)  /h  , q = 1 ~ P and 
f(t.),  f'Ct^),  f"(t^)  are  the  position,  velocity  and  acceleration  components 
at  time  t.  . 

l 

The  satellite's  velocity  (x  , y , z ) immediately  after  the  gross  attitude 

m m m 

manoeuvre  is  given  by 


where  (6vT#  dv^,  <5vR)  are  the  transverse,  normal  and  radial  components  of  delta- 
velocity  resulting  from  the  manoeuvre. 
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Svanmary 

Language 

Author 

Function  statement 
Input  argument 
NAME 


Output  function 
INFIND 


FUNCTION  INFIND 

- The  function  finds  the  location  of  a named  area  on  a 
specified  disc  file. 

- PLAN  for  use  with  1900  Fortran. 

- A. W.  Odell  (July  1973). 

- FUNCTION  INFIND  (NAME) . 


Name  of  an  area  on  the  disc  file  up  to  12  characters  in 
length;  must  be  array  element  or  text  (Hollerith) 
constant. 


Number  specifying  an  area  on  the  disc  file,  0 if  the 
name  is  not  found  in  the  index. 


Use  of  COMMON  - The  first  128  integer  locations  of  blank  common  are 

used  as  temporary  working  space. 

Source  deck  - 37  cards,  including  3 comment  cards. 

Local  storage  used  - 37  words  for  program,  7 words  for  data. 

Subordinate  subprograms  - The  subroutine  UTD4. 

Explanation  - The  subroutine  assunes  that  an  index  has  been  set  up  on 

the  disc  file  using  subroutine  INITD  and  that  information  has  been  put  in  the 
index  using  subroutine  ADDINF.  Associated  with  each  name  in  the  index  is  a 
number  specifying  an  area  on  the  disc  file.  If  this  name  is  not  found,  INFIND 
is  set  to  0;  otherwise  it  is  set  to  the  associated  number. 

Remark : A Fortran  version  of  this  subroutine  exists. 
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FUNCTION  INTFRC 

- The  function  computes  the  integral  and  fractional  parts 
of  a number. 


Language 


Author 


Function  statement 


- PLAN,  for  use  with  1900  Fortran. 

- A.W.  Odell  (February  1971). 

- FUNCTION  INTFRC  (X). 


Input  and  output  argunents  - 


Output  function 


Real  number,  which  is  truncated  to  its  fractional  part. 


INTFRC 


Integral  part  of  X. 


Use  of  COMMON 


- None. 


Source  deck 


Local  storage  used 


- 73  cards  including  3 comment  cards  (ICL  code). 

- 18  words  for  program. 


Subordinate  subprograms  - None. 

Explanation  - X is  split  into  its  mathematical  integral  part  and  its 

fractional  part.  For  example: 

X - - 3.4 
I - INTFRC  (X) 

would  result  in  I being  set  to  -4  and  X to  0.6. 

A Fortran  version  of  this  subroutine  is  also  available. 
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Language 


Authors 


SUBROUTINE  IPTOCO 

- The  subroutine  computes  the  geocentric  cartesian  com- 
ponents of  the  position  and  velocity  of  an  earth  satel- 
lite in  PROP  axes,  given  the  date  and  time  and  a set 

of  the  standard  parameters  used  to  define  orbit 
injection. 

- USA  Standard  Fortran  (USAS  X3.9  - 1966). 

- G.E.  Cook  and  R.  Clarke  (October  1972). 


Subroutine  statement  - SUBROUTINE  IPTOCO  (X,  Y,  Z,  XDOT,  YDOT,  ZDOT,  V,  CA, 

AZ,  R,  XLAT,  XLONG,  MJD,  TIME). 


Input  arguments 


Speed,  v . 

Climb  angle,  6 . 


Azimuth  of  velocity  vector,  Y (measured  E of  N) . 
Radial  distance,  r . 


XLONG 


Geocentric  latitude,  $ . 
Longitude,  X . 

Modified  Julian  Day  number. 


Output  arguments 


Time  (fraction  of  a day)  such  that  t is  given  by 
MJD  + TIME  . 


X,Y,Z 

XDOT,  YDOT,  ZDOT 


XLONG 


Use  of  COMMON 


Source  deck 


Local  storage  used 


Geocentric  coordinates  of  satellite,  x,  y,  z . 
Geocentric  components  of  velocity  vector,  x,  y,  z . 
Longitude  in  inertial  axes,  X*  . 

- None. 

- 22  cards  (ICL  code). 

- 8 real  variables. 


Subordinate  subprograms  - None. 


Explanation  - Units  of  length  and  time  are  arbitrary,  except  for  the 

variables  MJD  and  TIME;  angles  are  in  radians. 
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The  components  of  position  and  velocity  are  evaluated  relative  to  the 
following  coordinate  system:  the  origin  0 is  at  the  earth's  centre  and  Oz  points 
towards  the  north  pole;  Ox  lies  in  the  plane  of  the  true  equator  of  date,  but 
instead  of  pointing  towards  the  true  equinox  of  date  it  points  towards  a projec- 
tion of  the  mean  equinox  of  the  epoch  1950.0  (MJD  33281.9234);  Oy  completes  the 
right-handed  system  Oxyz. 

The  injection  conditions  are  specified  relative  to  an  earth-fixed  system 
with  longitude  measured  eastwards  from  the  Greenwich  meridian.  The  longitude  A' 
relative  to  the  inertial  system  defined  above  is  found  by  adding  § to  A , 
where  0 is  the  'modified  sidereal  angle'  given  by 

0 - 100.075542  + 360.985612288(t  - 33282.0)  , 

i.e.  the  origin  has  been  adjusted  slightly  from  1950.0.  (The  modified  sidereal 
angle  differs  from  sidereal  time  only  because  of  the  choice  of  a non-standard 
reference  direction.) 

The  geocentric  components  of  position  and  velocity  are  obtained  from  the 
following  equations: 

x ■ r cos  $ cos  A' 

y ■ r cos  <p  sin  A' 


z - r sin  $ 


x - v |sin  0 cos  $ cos  A'  - cos  0 (cos  ¥ sin  $ cos  A + sin  ¥ sin  A)} 
y ■ v jsin  0 cos  + sin  A'  ♦ cos  0 (sin  ¥ cos  A - cos  ¥ sin  $ sin  A)} 


{sin  0 sin  + ♦ cos  0 cos  ¥ cos  $( 
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SUBROUTINE  LOTMAT 


Summary  - Given  a symmetric  and  positive  definite  matrix,  the 

subroutine  performs  triangular  decomposition  into  a 
lower  triangular  matrix,  which  when  postmultiplied  by 
its  transpose  gives  the  original  matrix. 


Language 

Authors 

Subroutine  statement 

Input  arguments 

V(L,  L) 

L 


- USA  Standard  Fortran  (USAS  X3.9  - 1966). 
~ G.E.  Cook  and  R.  Clarke  (October  1972) 

- SUBROUTINE  LOTMAT  (V,  C,  L) . 

- Matrix  to  be  decomposed. 

~ Number  of  rows  and  columns  in  V. 


Output  arguments 
C(L,  L) 

Use  of  COMMON 

Source  deck 


Local  storage  used 
Subordinate  subprograms  - None. 


- The  decomposed  lower  triangular  matrix. 

- None. 

- 24  cards,  including  3 comment  cards  (ICL  code). 

- 3 integer  variables. 


Explanation  ~ If  V is  any  symmetric  and  positive  definite  matrix 

with  elements  v^  , there  exists  a unique  lower  triangular  matrix  C , such 

that  V ■ CC'  , where  C'  is  the  transpose  of  C . The  elements  c..  of  C 
. lj 

are  determined  recursively  from  the  following  equations: 


r 
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SUBROUTINE  MODAT 

r** 

Summary 

- The  subroutine  evaluates  upper-atmosphere  density, 
density  scale  height  and  scale  height  gradient  using 
a simple  analytic  model. 

Language 

- 1900  Fortran. 

Authors 

- G.E.  Cook  and  K.J.  Tomlinson  (April  1969). 

• 

Subroutine  statement 

- SUBROUTINE  MODAT  (HEIT,  TINF,  TGRADO,  DEN,  SCALHT,  SHGRAD) . 

Input  arguments 

• 

HE  IT 

Height  above  the  earth's  surface,  y . 

TINF 

Exospheric  temperature,  T 

00 

TGRADO 

Atmospheric  temperature  gradient  dT/dy  at  the 
reference  altitude,  y^  (120km). 

Output  arguments 

DEN 

Atmospheric  density,  p (g/cm  ). 

SCALHT 

Density  scale  height  (km). 

SHGRAD 

Density  scale  height  gradient. 

Use  of  COMMON 

- None. 

Source  deck 

- 46  cards,  including  2 comment  cards  (ICL  code). 

Local  storage  used 

- 24  real  array  elements,  16  real  variables,  1 integer 

variable. 

Subordinate  subprograms 

- None. 

; 

• 

Explanation  - The  values  of  density  and  densitv  scale  heieht  are 

obtained  from  a simple  analytic  model  of  the  Earth's  upper  atmosphere.  If  gQ 
denotes  the  local  value  of  the  acceleration  due  to  gravity,  the  geopotential 
height  above  the  reference  altitude  y^  is  defined  by 

c - J |g(y)/g(y0)}dy  t 0) 

y0 

• 

#0  that 
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c - (y  - y0)<R  ♦ yn>/<R  ♦ y)  . 


(2) 


R being  Che  mean  radius  of  Che  EarCh.  The  CemperaCure  of  Che  aCmosphere  is 
represenced  as  a funcCion  of  geopoCenCial  heighc  by  Che  expression 


T(y)  - Tjl  - a exp (-  tO\  , 


(3) 


where  is  Che  exospheric  CemperaCure  and  a and  t are  consCanCs  defined 

by 

T(y0> 


a - 1 - 


(4) 


and 


T - 


T - T(y)  \dy/ 
oo  70  7 y*yr 


(5) 


If  Che  aCmosphere  is  assumed  Co  be  in  diffusive  equilibrium  above  Che  reference 
altitude,  Che  number  c 
mass  nu  is  given  by 


alcicude,  Che  number  densicy  n^  of  Che  ich  consCiCuenC  of  molecular  (or  aComic) 


dn. 


n.  dy 
l 


mi8  1 dT  ( v 

kT  T dy  ( a) 


(6) 


where  k is  BolCzmann’s  constant  and  a is  the  thermal  diffusion  factor. 

Uich  Che  temperature  profile  (3),  equation  (6)  can  be  integrated  Co  give 


.l+Y.+o 


ni  " ni(y0>  jl  -V  exp(-  xc>{  eXp(_  YiTC>  ’ 


where  Y.  - 


mig(y0) 


i tkT 


k 


I 
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The  density  p is  given  by 


P 


z 


n.m. 
1 1 


For  the  constant  boundary  conditions  at  the  reference  altitude  of  120km 


we  use  the  values  assumed  by  Jacchia^ 
sion  profiles: 

T ( 1 20)  - 

n(N2)  ■ 

n(02)  - 

n(0)  - 

n(He)  - 


in  the  construction  of  his  static  diffu- 

355  K 

4.0  x io"  cm-3 

7.5  x 1010  cm-3 

7.6  x io10  cm-3 
3.4  x 10^  cm  3 


For  hydrogen  we  also  follow  Jacchia  and  take  the  concentration  at  500km  to  vary 
with  according  to  the  relation 

log1()  n(H;  500)  - 73.13  - 39.40  log,0  ♦ 5.5  (log10  TJ2  . 

The  thermal  diffusion  factor  a is  taken  as  -0.4  for  helium  and  as  zero  for  the 
other  constituents. 


k 
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SUBROUTINE  0RBIT2 

- The  subroutine  integrates  a satellite  orbit  to  a given 
epoch  and  interpolates  to  find  the  geocentric  position 
and  velocity  components.  If  required,  the  integration 
may  be  continued  to  a second  epoch  and  another  inter- 
polation performed.  The  effects  of  a gross  attitude 
manoeuvre  at  an  intermediate  apse  may  be  included. 

- 1900  Fortran. 


Author  - M.D.  Palmer  (August  1975) 

Subroutine  statement  - SUBROUTINE  ORBIT 2 


Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CINTEG/ , 

/CONST/,  /CSM00N/  and  /INTP/  are  used  as  follows:- 

Input  arguments  in  /CINTEG/  - 

PA(3)  Satellite's  geocentric  acceleration  components. 

* 2 

v ■ (x,  y,  z)  at  the  latest  integration  step  (km/s  ). 

Input  arguments  in  / CONST/  - 

3 ~2 

EMU  The  earth's  gravitational  constant  (km  s ). 

Input  arguments  in  /CSMOON / - 

MJDT  Modified  Julian  day  number  of  the  current  time. 

TIMET  The  current  time,  in  fractions  of  a day,  relative 

to  MJDT. 


Input  arguments  in  /INTP/  - 


EP 


TINT  I 


The  time  of  epoch,  in  fractions  of  a day,  relative 
to  MJDOCH. 

The  time  at  which  the  first  interpolation  is  required, 
measured  in  days  from  epoch. 


The  time  at  which  the  second  interpolation  is  required, 
measured  in  days  from  epoch.  If  TINT2  * 0.0,  the 
integration  is  terminated  after  first  interpolation. 

The  number  of  the  apse,  at  which  the  manoeuvre  is 
required.  NAM  ■ 0 if  no  manoeuvre  is  to  be  included. 
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Input  and  output  arguments  in  /CINTEG/  - 


MJDOCH 


PP(3) 


PV(3) 


On  input,  the  Modified  Julian  day  number  of  epoch. 
On  output  the  Modified  Julian  day  number  of  the 
apse  at  which  the  attitude  manoeuvre  was  performed. 

On  input  the  satellite's  geocentric  position  com- 
ponents, r.  ■ (x,  y,  z),  at  epoch.  On  output,  the 
satellite's  position  components  at  the  apse  at 
which  the  attitude  manoeuvre  was  performed  (km) . 

On  input  the  time,  in  seconds, of  epoch  relative 
to  MJDOCH.  On  output,  the  time,  in  seconds, of 
the  manoeuvre  apse  relative  to  MJDOCH. 

On  input,  the  satellite's  geocentric  velocity 
components,  v * (x,  y,  z)  at  epoch.  On  output, 
the  satellite's  velocity  components  immediately 
after  the  manoeuvre  (km/s). 


Output  arguments  in  /CINTEG/  - 


Source  deck 


Local  storage  used 


Subordinate  subprograms 


Explanation  - If  a gross  attitude  manoeuvre  is  to  be  included, 

subroutine  GATMAN  is  called  to  find  the  time  of  the  manoeuvre  and  the  satel- 
lite's geocentric  position  and  velocity  components  immediately  after  it  has 
taken  place. 

The  subroutine  starts  the  integration  of  the  orbit  either  from  the 
manoeuvre  time  or  from  the  initial  epoch  by  calling  subroutine  DEQRSPT  with  the 
statement  CALL  DEQRSPT  (IND,  PDAUXP) . On  tha  first  call  IND  ■ 1 and  on  the 
second  and  subsequent  calls  IND  ■ 2. 


The  initial  values  of  dx/ds,  dy/ds  and  dz/ds. 


The  initial  value  of  dt/ds  (s  ). 

- 82  cards  (ICL  code) . 

- 15  real  array  elements,  25  real  variables  and 
2 integer  variables. 

- The  subroutines  DEQRSPT,  GATMAN,  IPTOCO,  MODAT, 
PDAUXP,  SMPOS , TRINV  and  the  functions  ANGLE, 
ANGLI , INFIND,  INTFRC , SOLVIN,  SCPROD  and  UTD4 . 
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After  each  integration  step,  the  time  elapsed  from  epoch  (T2)  is  calcu- 
lated and  compared  with  TINT1.  If  T2  < TINT1  the  time,  position,  velocity  and 
acceleration  components  are  stored  and  another  integration  step  performed. 

(To  reduce  computation  time,  this  storage  is  not  carried  out  if  T2  < 0.9  TINT 1 . ) 

When  T2  i TINT1,  an  interpolation  is  carried  out  using  the  equations:- 
f(t)  - q3|jl  + 3p  + 6p2)  f(t  j)  + hp(l  + BpJf'Uj)  + | p2f" (tj)J 

+ P3[(1  + 3q  + 6q2)f(t2)~  hq(l  + 3q)f  ' ( t^  + ^ q2f » (t^] 
for  position  components  and 

f'(t)  = 30h_1p2q2[f(t2)  - f(tj)]  + q2( 1 + 5p) ( 1 - 3p)f ' (t j ) 

+ p2(l  + 5q) ( I - 3q)f ' (t2)  + | pq2 (2  - 5p)f "<t j) 

- | p2q(2  - 5q)f"(t2) 

for  velocity  components,  where  h ■ t2  - tj  , P“(t-t)/h,q*l-p  and 

f(t^),  f'(t^),  f'(t^)  are  the  position,  velocity  and  acceleration  components 

at  time  t.  . 

l 

f (t)  and  f * (t)  are  then  written  to  an  unformatted  disc  file  on  channel  2. 
If  TINT2  ^ 0.0,  the  above  process  is  repeated. 
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Language 

Author 

Subroutine  statement 


SUBROUTINE  PDAUXP 

- The  subroutine  calculates  the  cartesian  components 

of  the  geocentric  accelerations  due  to  the  gravitational 
attractions  of  the  earth,  sun  and  moon,  atmospheric 
drag  and  the  precession  of  the  earth's  polar  axis.  It 
will  also,  if  required,  evaluate  the  partial  deriva- 
tives of  the  acceleration  components  with  respect  to  the 
initial  components  of  the  position  and  velocity  at 
epoch. 

- 1900  Fortran. 

- M.D.  Palmer  (November  1974). 

- Subroutine  PDAUXP. 


Use  of  COMMON  - Certain  quantities  in  the  comnon  blocks  /CINTEG/ , /CON/, 

/CONST/,  /CSMOON/  and  /PETURB/  are  used  as  follows: 

Input  arguments  in  /CINTEG/  - 

MJDOCH  Modified  Julian  day  number  of  epoch. 

NEQ  Number  of  equations  being  integrated  (22  if  partial 

derivatives  of  acceleration  are  required,  4 otherwise). 


IAUX 
X,  Y,  Z 
PD ( 1 8) 

XVEL,  YVEL,  ZVEL 
T 

PDVEL( 1 8) 


Flag,  see  explanation. 

Cartesian  position  components  (x,y,z)  - r . 

Partial  derivatives  of  position  with  respect  to  the 
initial  position  (r^)  and  velocity  (v^) ; 3r/3(r0,VQ)  . 

The  latest  values  of  dx/ds,  dy/ds  and  dz/ds  . 

Time,  t , relative  to  MJDOCH  (seconds). 

Derivatives  with  respect  to  the  independent  variable 

of  position 


8 of  the  partial  derivatives 


_d_ 

ds 


3r 


3(r0*V 


Appendix  B 

Input  arguments  in  /CON/  - 
EJ5,  EJ6,  EJ7,  EJ8,  EJ9  Coefficients 

>J j , • • • » 

Fully  normalized  coefficients  of  certain  tesseral 
harmonics,  viz.  C^,  S^,  C ~,  S^\  C^,  . 

Input  arguments  in  /CONST/  - 

Earth's  gravitational  constant,  . 

Coefficients  of  the  earth's  zonal  harmonics,  J2,  J^, 
J4  • 

Fully  normalized  coefficients  of  certain  tesseral 
harmonics,  viz.  C^,  S^,  cJJ\  S^\  Cj^,  S^,  C^, 

S42»  C44»  S44  * 

Mean  equatorial  radius  of  the  earth,  R . 

Mean  rotation  rate  of  the  earth's  polar  axis,  in  rad/s. 
Sun's  gravitational  constant,  . 

Moon's  gravitational  constant,  . 

Input  arguments  in  /CSMOON/  - 

XS,  YS,  ZS  Cartesian  components  of  the  sun's  position  (r  ) at 

s 

the  current  time,  computed  if  IAUX  < 0 . 

XM,  YM,  ZM  Cartesian  components  of  the  moon's  position  (r  ) at 

m 

the  current  time,  computed  if  IAUX  < 0 . 

Input  arguments  in  /PETURB/  - 


SMEAN 

Mean  value  over  three  days  of  the  solar  10.7cm  radia- 
tion flux  in  units  of  10  m ^ Hz  * . 

SOLBAR 

Mean  value  over  three  solar  rotations  of  the  10.7cm 

-22  -2  -1 

radiation  flux  in  units  of  10  W m Hz  . 

RATE 

Angular  velocity  of  the  atmosphere,  in  rad/s. 

ARMA 

The  product  AM  x HALF CD  x 1000  where  AM  is  the  satel- 
lite's area-to-mass  ratio  and  HALFCD  is  the  quantity 
JCD  where  CD  is  the  satellite's  drag  coefficient. 

If  ARMA  - 0.0  , drag  terms  are  excluded. 


EMU 


EJ2,  ej^ 


C22,  S22 
C31 , S31 
C33,  S33 
C42,  S42 
C44,  S44 

ERAD 


EOMEGA 

SUNMU 

SELMU 


C32,  S32] 
C41 , S41  > 
C43,  S43  J 


of  the  earth's  zonal  harmonics, 
J9  * 
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TNITE  Minimum  night-time  temperature  (K) . 

LSP  If  LSP  > 0 , luni-solar  perturbations  are  included. 

NSR  If  NSR  > 0 , solar  radiation  pressure  terms  are 

included  but  see  explanation. 

Output  arguments  in  /CINTEG/  - 

XACC,  YACC,  ZACC  Values  of  d2x/ds2,  d2y/ds2,  d2z/ds2  at  the  current 

time. 


PDACC( 1 8) 


Value  of  dt/ds  at  the  current  time. 

Second  derivatives  with  respect  to  the  independent 
variable  s of  the  partial  derivatives  of  position 


31 

^ 2 3(rn,v.) 


XVELT , YVELT,  ZVELT  Cartesian  components  of  velocity,  x,  y,  z,  at  the 


current  time. 


XACCT,  YACCT,  ZACCT  Cartesian  components  of  acceleration,  x,  y,  z,  at  the 


current  time. 


PDACCT ( 1 8) 


Partial  derivatives  of  acceleration  with  respect  to 
the  initial  position  and  velocity  at  epoch 


dt2  3(r0‘V 


Output  arguments  in  /CSMOON/  - 


TIMET 


Source  deck 


Modified  Julian  day  number  of  the  current  time. 

The  current  time,  in  fractions  of  a day,  relative  to 
MJDT . 

- 155  cards  (ICL  code). 


Local  storage  used  - 51  real  variables  and  3 integer  variables. 

Subordinate  subprogram  - The  subroutines  MODAT,  SMPOS  and  TRINV  and  the  functions 

ANGLE,  ANGLI,  INFIND  and  SCPROD. 

Explanation  - When  the  subroutine  is  called  for  the  first  time,  or  for 

a change  of  epoch,  IAUX  must  be  set  to  -1.  For  subsequent  entries,  when  the 
time  has  changed  since  the  previous  entry,  IAUX  should  be  set  to  zero.  If  the 


1 

Appendix  B — 
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time  has  not  changed  IAUX  may  be  set  to  1 to  prevent  recomputation  of  certain 
terms.  (All  these  cases  are  allowed  for  by  the  calling  segment  in  POINT.) 

The  cartesian  components  of  the  geocentric  acceleration  acting  on  the 
satellite  are  found  by  adding  contributions  from  the  gravitational  attractions 
of  the  earth,  sun  and  moon,  atmospheric  drag  and  the  precession  of  the  earth's 
polar  axis.  The  subroutine  is  structured  so  as  to  permit  the  easy  insertion  of 
solar  radiation  pressure  at  a future  date. 

The  earth  harmonic  terms  are: 


I H « . j£>y 

n,m  ' L - 


where  the  sidereal  angle  0 - 0q  + , 

W - (x  + jy)/re^0)  , 

(m ) . 

is  the  mth  derivative  of  the  Legendre  polynomial,  P^(z/r) 

E0,0  “ 1 ' 

En,0  " ' Jn 


and 


E «■  [2(2n  + 1 ) (n  - m)l/(n  ♦ m)ll^(C  - jS  ) 
n,m  • n,ra  n,m 

Terms  are  included  up  to  the  ninth  zonal  harmonic  and  the  4,4  tesseral 
harmonic. 

The  accelerations  due  to  the  gravitational  attractions  of  the  sun  and 
moon  are  given  by 


<rb  -ll3  - rb/lrb' 


b being  's'  for  the  sun  and  'm'  for  the  moon. 

Atmospheric  drag  is  included  if  ARMA  =£  0 . A modified  version  of  the 
Jacchia  1965  model7  is  used  to  find  the  ambient  air  density.  Firstly  the  sub- 
routine computes  the  satellite's  height  and  latitude  ($),  the  declination  of  the 
sun  (6)  and  the  hour  angle  (H)  of  the  satellite  relative  to  the  sun.  The 
exospheric  temperature  (T  ) is  determined  using  the  equation, 
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T - T, 


NITE 


1 + 0.280  ♦ 0.28  [cos  {-^}2‘5  - e]  [cos  , 


where  0 - sin  ^ 


2.5 


and 


- [h  - 0.78539816  + 0.20943951  sin  [l  + 0.785398I6]J  . 


The  atmospheric  temperature  gradient  at  the  reference  altitude  (120km)  is 
given  by: 

Tgrad0  " (T»  ’ 355) (0.029  exp[-  x2/2]) 


T - 800 


where  x ** 


750  + 1.722  x 10-4(T  - 800)2 


Subroutine  MODAT  is  called  to  determine  the  atmospheric  density  p . 

The  force  acting  on  the  satellite  is  given  by  Jp|v|2SCD  where  CQ  is  the  drag 
coefficient,  S is  the  effective  cross-secticnal  area  perpendicular  to  the  air 
flow  and  V is  the  velocity  of  the  satellite  relative  to  the  ambient  air.  V is 
given  by 

V « v - ojg  x £ 

where  £ and  v are  the  position  and  velocity  vectors  of  the  satellite. 

If  the  cartesian  components  of  V are  V , Vv,  V then  the  contributions 
to  the  acceleration  components  are: 


x 


etc. 


where  M is  the  mass  of  the  satellite. 

SCD 

The  variable  ARMA  is  the  quantity  x 1000  , the  1000  being  a unit 

conversion  factor. 


The  contribution  to  the  acceleration  due  to  the  precession  term  is 


p(xz  - zx)  . 
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Summary 

Language 

Author 

Function  statement 
Input  arguments 
A(  I ) , B (1 ) 

NA,  NB 

N 

Output  function 
SCPROD 


Use  of  COMMON 
Source  deck 


FUNCTION  SCPROD 

- The  function  gives  the  scalar  (inner)  product  of  two 
arrays . 

- PLAN,  for  use  with  1900  Fortran. 

- A.W.  Odell  (March  1973). 

- FUNCTION  SCPROD  (A,  B,  NA,  NB,  N) . 


Locations  of  the  first  elements  to  be  multiplied  (must 
be  array  elements). 

Increments  in  the  arrays.  A,  B of  elements  to  be 
multiplied. 

Number  of  elements  to  be  multiplied. 


The  scalar  product 


N 

^ A ( 1 + (i  - I )NA) B( 1 + (i  - I )NB) . 
i-1 


- None. 


- 25  cards,  including  1 comment  card  (ICL  code). 


i 


Local  storage  used  - 2 words  for  variables,  plus  22  words  for  program. 

Subordinate  subprograms  - None. 


Explanation  - The  arrays  in  the  calling  routine  may  have  any  dimen- 

sions although  they  are  treated  as  one-dimensional  arrays  in  the  function.  For 
simplicity  A and  B are  dimensioned  * 1 * . This  may  cause  the  function  to  fail 
in  TRACE 2,  so  it  should  not  be  compiled  in  this  mode.  If  N < 0 , the  function 
will  return  zero. 


j 
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Summary 

Language 

Author 

Subroutine  statement 
Dummy  arguments 


Appendix  B 


SUBROUTINE  SMPOS 

- The  subroutine  computes  the  geocentric  cartesian 
coordinates  (km)  of  the  sun  and  moon  at  a given  time. 

- ASA  Fortran  (Standard  Fortran  4). 

- A.W.  Odell  (Hay  1973). 

- SUBROUTINE  SMPOS. 

- None. 


! 

i 


Use  of  COMMON  - Certain  quantities  in  common  block  /CSMOON/  are  used 

as  follows: 

Input  arguments  in  /CSMOON/  - 

MJDT  Modified  Julian  day  number  of  the  required  time. 

TIMET  Time,  as  a fraction  of  a day,  since  0 hours  on  MJDT. 

Output  arguments  in  /CSMOON/  - 

XS,  YS,  ZS  The  cartesian  coordinates  (x  ,y  ,z  ) = r of  the  sun. 

s’^s’  s s 

XM,  YM,  ZM  The  cartesian  coordinates  (x  ,y  ,z  ) = r of  the  moon. 

ra’-'m'  m m 

Input  and  output  argument  in  /CSMOON/  - 
TABLE (43)  Working  space. 

Source  deck  - 27  cards  including  1 comment  card  (ICL  code) . 

Local  storage  used  - 2 integer  variables,  1 logical  variable. 

Subordinate  subprograms  - The  subroutine  UTD4,  and  the  functions  INFIND  and 

SCPROD. 


Explanation  - The  subroutine  obtains  the  sun-moon  coordinates  by 

interpolation  in  a table  of  daily  positions,  with  2nd  and  4th  differences,  using 
Everett's  interpolation  formula.  The  table  is  read  from  a disc  file  when 
required,  using  the  subroutine  UTD4.  It  is  stored  in  an  area  called  SUNMOONTABLE , 
the  location  of  which  is  found  using  the  function  INFIND.  The  area  is  contained 
in  a disc  file,  which  must  have  been  opened,  before  entry  to  SMPOS,  on 
channel  7. 


The  table  is  stored  as  follows:  firstly,  the  modified  Julian  day  numbers 

of  the  first  and  last  sets  of  coordinates  (2  integers),  then  data  for  each 

2 2 A A 

midnight  as  follows  (sets  of  18  reals):  r,r,Ar,Ar,Ar,  A r . 

6 s*  m*  s’  m’  s’  m 


i 


J 
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if  fQ»  f|  denote  values  of  f at  t « 0 and  t ■ 1 , Everett's  inter- 
polation formula  to  the  4th  differences11  gives: 

f(t)  - DQ  + t^D,  + (t  - 1)(d2  + (t  + 1)(d3  + (t  - 2)(d4  + (t  + 2)D5))j)  , 

where  - 62nfQ/(2n)!,  D2n+,  - 62n(f,  - fQ)/(2n  +1)1. 

The  data  was  originally  obtained  on  punched  cards  from  JPL12,13and  has  been 
transformed,  before  storing  on  the  disc,  into  the  SAO/PROP14  system  of  axes, 
using  the  subroutine  AX1950.  If  the  time  input  to  SMPOS  lies  outside  the  range 
of  data  stored  on  the  disc,  a STOP77  statement  will  be  obeyed. 
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FUNCTION  SOLVIN 


Summary 


Language 

Author 

Function  statement 
Input  arguments 


- Given  a function  and  its  derivative  at  two  points,  the 
function  solves  four  problems  involving  cubic 
interpolation. 

- ASA  Fortran  (Standard  Fortran  4). 

- A.W.  Odell  (January  1975). 

- FUNCTION  SOLVIN  (MODE,  M,  FO,  FI,  FDO,  FD1 , F) . 


MODE 

H 

FO,  FI 
FDO,  FD1 

F 

Output  function 


Number  specifying  problem  to  be  solved,  see  explanation. 

Difference,  h , between  the  arguments  of  the  function 
(i.e.  h - x(  - xQ). 

Function  values  f^  at  x - x^  and  fj  at  x ■ x^  . 

Derivative  values  f^  , at  x ■ x^  and  fj  at 
x - x(  . 

Required  function  value  or  argument  - see  explanation. 


SOLVIN 

Use  of  COMMON 
Source  deck 
Local  storage  used 


Solution  to  problem  - see  explanation. 

- None. 

- 33  cards,  including  5 comment  cards  (ICL  code). 

- 9 real  variables. 


Subordinate  subprograms  - None. 

2 3 

Explanation  - A cubic  polynomial  P(x)  ■ f^  ♦ f^x  + + a^x  is 

first  fitted  to  the  data  giving  ■ (3B  - A)/h  and  a^  - (A  - 2B)/h2  where 
A - f*  - fQ  and  B - (f,  - fQ)/h  - f^  . 

Then, 

(1)  if  MODE  - 1 , Newton's  method  is  used  to  find  the  value  of  x for  which 

P(x)  - F , 

(2)  if  MODE  ” 2 , Newton's  method  is  used  to  find  the  value  of  x for  which 

P'(x)  - F , 

(3)  if  MODE  - 3 , P(F)  is  evaluated,  and 

(4)  if  MODE  - 4 , P'(F)  is  evaluated. 
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Summary 

Language 

Author 

Subroutine  statement 
Input  arguments 
Y 
X 

Output  arguments 
R 

TH 

Use  of  COMMON 
Source  deck 
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SUBROUTINE  TRINV 

- The  subroutine  obtains  polar  coordinates  from  (two- 
dimensional)  cartesian  coordinates. 

- ASA  Fortran  (Standard  Fortran  4). 

- R.H.  Gooding  (May  1968). 

- SUBROUTINE  TRINV  (Y,  X,  R,  TH) . 

Cartesian  y-coordinate  (arbitrary). 

Cartesian  x-coordinate  (arbitrary). 

Polar  r-coordinate. 

Polar  9-coordinate. 

- None. 

- 8 cards,  including  2 comment  cards  (ICL  code). 

- None. 


Local  storage  used 
Subordinate  subprograms  - None. 

Explanation  - The  equations  r cos  0 = x , and  r sin  0 * y are 

solved  for  r and  0 , using  the  standard  ATAN2  function  to  give  0 between 
-it  and  +n  . If  x « y - 0 the  ATAN2  function  is  not  used  and  0 is  set  to 
zero. 

Remarks 

(1)  The  subroutine,  if  used  twice,  provides  a convenient  solution  of  the  three- 
dimensional  cartesian-to-polar  transformation.  Thus  to  solve  the  equations, 

r sin  9 cos  $ - x , 

r sin  0 sin  $ - y , 

and 

r cos  0 - z 

for  r,  0 and  $ , the  following  two  statements  will  suffice: 

CALL  TRINV  (Y,  X,  RSINTH,  PHI)  , giving  r sin  0 and  * 

CALL  TRINV  (RSINTH,  Z,  R,  TH)  , giving  r and  0 . 


and 


64 

Moreover,  this  solution  will  give  maximum  accuracy  for  0 and  $ 
angles  set  in  the  correct  quadrant. 

(2)  The  actual  input  and  output  arguments  must  be  distinct. 
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ANGLE 


ANGL1 


(a) 

(a) 


(b) 


ANST 
APMAN 
ARAN0R2 


(c) 

(c) 


ARAN0R6 
BLOCK  DATA 
COTOEL 
DR(d) 

DEQRSPT 
DVST(b) 


(a) 


EAFKEP 

FPA(<1) 


INFIND 

INTFRC 

INTPTB 

LOTMAT 

MATADD 


(a,e) 

(a,e) 

(a) 


MATMUL 
.(a) 


MODAT 


PDAUXP 

OFPA(b) 


(a,f  ) 


RANDOMNO(g) 

RLEASE(g) 

SCPROD(a,e) 

SIDANG 

SMPOS(a) 

STATIS 

TRINV(a) 

TRNFRM 

XROT 

UTD4(g) 


Appendix  C 

SYNMAP  PROGRAM  UNITS 

reduces  an  angle  to  the  range  0 to  2tt 

reduces  an  angle  to  the  range  -it  to  n 

performs  station  acquisition  by  the  shortest  path 

determines  the  ABM  burn  point  and  the  nominal  ABM  firing  direction 

generates  a pair  of  random  normal  deviates 

generates  six  random  normal  deviates 

sets  certain  constants 

converts  coordinates  to  osculating  elements 

Finds  the  nominal  ABM  firing  direction  for  a specified  drift 
orbit  drift  rate 

integrates  second  order  differential  equations 

performs  station  acquisition  with  minimum  delta-velocity 
expenditure 

determines  eccentric  anomaly  from  Kepler's  equation 

finds  the  nominal  ABM  firing  direction  for  a specified  flight 
path  angle 

finds  a named  area  on  a disc  file 

converts  a number  into  its  integer  and  fractional  parts 

computes  the  position  and  velocity  components  by  interpolation 

performs  triangular  decomposition  of  a matrix 

performs  matrix  addition 

performs  matrix  multiplication 

computes  atmospheric  density 

auxiliary  for  DEQRSPT,  computes  perturbing  accelerations 

optimises  the  flight  path  angle  and  finds  the  corresponding  nominal 
ABM  firing  direction 

generates  a random  number  in  the  range  0 to  1 . 
releases  peripheral  unit 
forms  scalar  product 

calculates  the  sidereal  angle  at  a given  epoch 
reads  sun/moon  coordinates  from  disc  and  interpolates 
produces  statistical  information 

determines  polar  coordinates  from  cartesians  (two-dimensional) 
performs  the  transformation  X ■ A + B C 
rotates  rectangular  axes  through  a given  angle 
permits  disc  core  transfers. 


1 
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(a)  These  subprograms  are  common  to  both  P0INT2  and  SYNMAP.  Specifications 
for  them  are  given  in  Appendix  B. 

(b)  These  subroutines  appear  in  the  subprogram  that  calls  them  as  STAC. 

(c)  These  subroutines  are  copies  of  ARANOR  for  which  a specification  is  given 
in  Appendix  B.  Only  the  names  have  been  changed.  Groups  of  six  random  normal 
deviates  are  only  required  for  runs  which  take  tracking  errors  into  account. 

Separate  subroutines  are  therefore  provided  so  that  the  same  pairs  of  random 
normal  deviates  will  be  generated  irrespective  of  whether  tracking  errors  are  to 
be  included  or  not. 

(d)  These  subroutines  appear  in  the  subprogram  that  calls  them  as  DUM. 

(e)  These  subprograms!  written  in  PLAN , are  provided  as  semicompiled  segments. 

(f)  PDAUXP  appears  in  the  subprogram  that  calls  it  as  DAUX. 

(g)  Provided  automatically  by  1900  series  compilers  and  not  described  here. 
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SUBROUTINE  ANSI 

Summary  - The  subroutine  calculates  the  delta-velocity  required 

for  a geosynchronous  satellite  to  acquire  station  by 
the  shortest  path  after  ABM  burn. 

Language  - 1900  Fortran 

Author  - M.D.  Palmer  (August  1975) 

Subroutine  statement  - SUBROUTINE  ANST(NO) 

Input  argument 

NO  The  case  number.  If  NO  a 0,  no  line  printer  output  is 

produced . 

Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CNST/, 

/STATION/  and  /STATS/  are  used  as  follows :- 

Input  argument  in  /CNST/  - 

ENO  Earth's  rotation  rate  (rad/d). 

RAD  Conversion  factor,  degrees  to  radians. 

Input  arguments  in  /STATION/  - 

X,  Y,  Z Satellite's  geocentric  position  components, 

r^  ■ (x,  y,  z)  at  ABM  burn  (km). 

XDOT,  YDOT,  ZDOT  Satellite's  geocentric  velocity  components  at 

yd  * (x,  y,  z)  immediate iy  after  ABM  burn  (km  s 

Modified  Julian  day  number  corresponding  to  the  ABM 
burn  time. 

Time  of  ABM  burn,  in  fractions  of  a day,  relative  to 
MJDB. 

The  time  allowed  for  tracking  and  orbit  determination 
after  ABM  burn  (days). 

Cl  ♦ TR  is  the  total  time  available  for  station  acquisi- 
tion (days). 

Station  longitude  (rad  E of  Greenwich). 

Minimum  and  maximum  permitted  values  of  drift  orbit 
inclination  (rad). 


MJDB 


TB 


TR 


I 


Cl 


PHI 


XIL,  XIH 


l 
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Input  arguments  in  /STATS/  - 


RAS,  DEC 


Right  ascension  and  declination  of  the  actual  ABM 
firing  direction  (rad) 


Output  arguments  in  / STATS/  - 


XINCD 


XLOND 


Source  deck 


Local  storage  used 


Total  delta-velocity  required  for  station  acquisition 
(ms  1 ) . 

Drift  orbit  inclination  (degree). 

Right  ascension  of  the  drift  orbit  ascending  node 
(degree) . 

Satellite's  longitude  at  the  ABM  burn  point  (°E). 
Drift  orbit  semi  major  axis  (km). 

Drift  orbit  eccencricity. 

- 50  cards  (ICL  code) . 

- 2 real  array  elements;  21  real  variables,  and 
1 integer  variable. 


Subordinate  subprograms  - The  subroutines  COTOEL  and  TRINV  and  the  functions 


ANGLE,  EAFKEP  and  SIDANG. 


Exp lanation  - The  subroutine  uses  linearised  equations  to  calculate 

the  total  delta-velocity  required  for  a synchronous  orbit  satellite  to  acquire 
station  by  the  shortest  path  from  the  ABM  burn  point.  The  total  delta-velocity 
required  for  this  strategy  is  not  necessarily  a minimum. 

The  value  of  the  earth's  gravitational  constant  is  taken  to  be 
3 2 

398616.82  km  /s  which  takes  into  account  the  effect  of  the  second  zonal 
harmonic  Jg  • 

COTOEL  is  called  to  find  the  osculating  elements  of  the  drift  orbit 
immediately  after  ABM  burn.  The  inclination  is  compared  with  the  maximum  and 
minimum  permitted  values,  and  if  necessary,  the  delta-velocity  required  to 
place  the  inclination  within  these  bounds  is  calculated.  The  direction  in 
which  the  satellite  is  to  be  moved  is  identified  and  the  necessary  drift  rate 
computed.  If  the  actual  drift  rate  is  inadequate  or  in  the  wrong  direction,  a 
manoeuvre  is  performed.  I*  is  assumed  that  this  and  the  final  stopping 
manoeuvre  are  also  used  to  reduce  the  orbital  eccentricity.  Finally,  the  total 
velocity  increment  needed  to  remove  any  residual  eccentricity  is  calculated. 


T 
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If  NO  # 0,  two  lines  of  information  comprising  the  case  number,  the  ABM 
firing  attitude,  drift  orbit  elements,  longitude  of  the  ABM  burn  point  and  the 
required  delta-velocity  components  are  output  to  the  line  printer. 
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Summary 

Language 

Author 

Subroutine  statement 
Input  arguments 
DUM 

EP 

EPS 

GAM 

Use  of  COMMON 

Input  arguments  in  /ABM/ 
HX,  HY,  HZ 

DELTAV 

DRA 

SDIN 


SUBROUTINE  APMAN 

The  subroutine  simulates  the  firing  of  the  apogee 
boost  motor  (ABM)  during  the  launch  phase  of  a syn- 
chronous orbit  mission. 

1900  Fortran. 

M.D.  Palmer  (April  1975). 

SUBROUTINE  APMAN  (DUM,  EP,  EPS,  GAM). 


The  name  of  a subordinate  subroutine  which  determines 
the  ABM  firing  direction.  Depending  on  the  chosen 
ABM  firing  strategy  it  may  be:- 

FPA  for  a preset  flight  path  angle, 

OFPA  for  an  optimised  flight  path  angle, 

DR  for  a preset  drift  rate. 

Time  of  the  epoch  to  which  the  random  orbits  relate, 
in  fractions  of  a day,  relative  to  MJDOCH  (days). 

Epoch  time,  in  seconds  relative  to  MJDOCH 
(ie  EPS  - 86400  * EP) . 

If  DUM  = FPA,  GAM  is  the  required  flight  path  angle 
(rad) . 

If  DUM  e OFPA,  GAM  is  set  to  1000.0  and  used  as  a 
flag. 

If  DUM  ■ DR,  GAM  is  the  required  drift  rate  (rad  s '). 

Certain  quantities  in  the  common  blocks  /ABM/, 
/CINTEG/,  /CNST/ , /CONST/,  /CSMOON/  and  /STATION/ 
are  used  as  follows:- 


Direction  cosines  of  the  unit  vector  normal  to  the 
required  drift  orbit  plane. 

Nominal  ABM  velocity  increment  (km  s *) 

Desired  right  ascension  of  the  drift  orbit  (rad) 

sin  i,  , where  i.  is  the  desired  drift  orbit 
a a 

inclination. 
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Input  arguments  in  /CINTEG/  - 


MJDOCH 

P(3) 

V(3) 

A(3) 

Input  arguments  in  /CNST/  - 

RAD  Conversion  factor,  degrees  to  radians. 

Input  arguments  in  /CONST/  - 

3 -2 

EMU  Earth's  gravitational  constant  (km  s ). 

Input  arguments  in  /CSMOON/  - 

MJDT  Modified  Julian  day  number  of  the  current  time. 

TIMET  The  current  time  in  fractions  of  a day  relative  to 

MJDT. 


Modified  Julian  day  number  of  the  epoch  to  which  the 
random  orbits  relate. 

Satellite's  geocentric  position  components, 

r.  “ (x,  y,  z),  at  the  latest  integration  step  (km). 

Satellite's  geocentric  velocity  components, 
v ■ (x,  y,  z) , at  the  latest  integration  step  (km  s *). 

Satellite's  geocentric  acceleration  components, 

# # —2 

v = (x,  y,  z),  at  the  latest  integration  step  (km  s ). 


Input  arguments  in  /STATION/  - 

SUB  The  name  of  the  subordinate  subroutine  to  be  used  to 

determine  the  velocity  increment  required  for  station 
acquisition.  (Only  required  if  DUM  s OFPA.) 

Output  arguments  in  /ABM/  - 

PT(3)  Satellite's  geocentric  position  components, 

r^  ■ (x,  y,  *),  at  the  ABM  burn  point  (km). 

VT(3)  Satellite's  geocentric  velocity  components, 

v^  ■ (x,  y,  z) , at  the  ABM  burn  point  (km/s). 

RB  Satellite's  radial  distance,  r^,  at  the  ABM  burn 

point  (km) . 

VB  Satellite's  transfer  orbit  velocity,  v.  , at  *-he 

D 

ABM  burn  point  (km/s). 


i 


179 


Appendix  D 


73 


Output  argument!  in  /CINTEG/  - 


XVEL 

YVEL 

ZVEL 

TVEL 


The  initial  values  of  dx/ds,  dy/ds  and  dz/ds. 


The  initial  value  of  dt/ds 


). 


Output  arguments  in  /STATION/  - 

MJDB  Modified  Julian  day  number  corresponding  to  the  ABM 

burn  time. 

TB  Time  of  the  ABM  burn,  in  fractions  of  a day,  relative 

to  MJDB. 


Source  deck  - 66  cards  (ICL  code). 

Local  storage  used  - 8 real  variables;  9 real  array  elements  and  1 integer 

variable. 


Subordinate  subprograms  - The  subroutines  ANST,  COTOEL,  DEQRSPT , DR,  DVST,  FPA, 

INTPTB , MODAT,  OFPA,  PDAUXP,  SMPOS,  TRINV  and  the 
functions  ANGLE,  ANGL1 , EAFKEP,  INFIND,  INTFRC , SCPROD 
and  SIDANG. 

Explanation  - The  subroutine  is  provided  with  the  satellite's  geo- 

centric position  and  velocity  components  at  a point  on  the  transfer  orbit  some 
2i  hours  before  the  apogee  selected  for  ABM  burn.  The  orbit  is  integrated 
forward  to  the  point  at  which  the  transfer  and  drift  orbit  planes  intersect. 

At  this  point,  the  value  of  r_  . li  , where  r_  is  the  radius  vector  and  h is 
the  unit  vector  normal  to  the  required  drift  orbit  plane,  changes  sign. 

Thus  the  integration  steps  before  and  after  the  intersection  can  be  identified. 
The  time  of  intersection  (te  the  ABM  firing  time)  is  found  by  interpolation. 

A check  is  made  that  the  ABM  velocity  increment  (6v)  is  adequate, 

ie  v.  . h < flv  . If  not,  an  error  message  is  printed  and  a STOP  VB.H  instruc- 

— D — 

tion  obeyed. 

Depending  on  the  firing  strategy  specified,  a call  is  made  to  one  of  the 
subroutines  DR,  FPA  or  OFPA  to  find  the  nominal  ABM  pointing  direction.  The 
name  of  the  subroutine  to  be  used  appears  as  an  argument  in  the  call  statement. 
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BLOCK  DATA  (SYNMAP) 

Suninary  - The  Fortran  BLOCK  DATA  segment  is  used  to  set  initial 

values  for  certain  quantities  stored  in  common  blocks 
/CIOTEG/,  /CNST/,  /CON/,  /CONST/  and  /PETURB/. 

Language  ~ 1900  Fortran. 

Author  - M.D.  Palmer  (April  1975). 


Data  in  /CINTEG/ 


Variable  name 

Value 

1(2) 

4 

H 

it/48 

Data 

in  /CNST/ 

Variable  name 

Value 

PI 

3. 1415926536 

PI2 

1.5707963268 

ENO 

6.300387476 

RAD 

0.0174532925 

SR 

42164.75 

THEL 

65.0  l 

• 

THEH 

1 15.0  J 

Data 

in  /CON/ 

Variable  name 

Value 

EJ5 

-0.246E-6 

EJ6 

0.558E-6 

EJ7 

-0.326E-6 

EJ8 

-0-209E-6 

EJ9 

-0.094E-6 

C43 

1 .0390E- 

6 

S43 

-0.1192E- 

6 

C4I 

-0.5175E- 

6 

S41 

-0.4814E- 

6 

C32 

0.7783E- 

6 

S32 

-0.7552E- 

6 

Explanation 

The  number  of  equations  to 
be  integrated 

Integration  step  length 

Explanation 

IT 

tt/2 

Earth's  rotation  rate 
(rad/day) 

n / 180 

Synchronous  orbit  radius 
(km) 

Lowest  and  highest  permitted  values 
of  solar  aspect  angle  (degree) 

Explanation 

Earth's  fifth  zonal  harmonic  J,. 
Earth's  sixth  zonal  harmonic  J, 

D 

Earth's  seventh  zonal  harmonic 

Earth's  eighth  zonal  harmonic  Jg 

Earth's  ninth  zonal  harmonic 

Tesseral  harmonic  coefficient  C. , 

43 

Tesseral  harmonic  coefficient  S,_ 

43 

Tesseral  harmonic  coefficient  C. , 

41 

Tesseral  harmonic  coefficient  S, , 

41 

Tesseral  harmonic  coefficient  C^ 
Tesseral  harmonic  coefficient  S^ 
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Data  in  /CONST/  - 

Variable  name 

Value 

Explanation 

EMU 

398601.3  km3  s“2 

Earth's  gravitational  constant 

EJ2 

1082.637E-6 

Earth's  second  zonal  harmonic  ^ 

EJ3 

-2.5310E-6 

Earth's  third  zonal  harmonic 

EJ4 

-1.6190E-6 

Earth's  fourth  zonal  harmonic 

C22 

2.4369E-6 

Tesseral  harmonic  coefficient  C ^ 

S22 

-1.4005E-6 

Tesseral  harmonic  coefficient  S^ 

C33 

0.7387E-6 

Tesseral  harmonic  coefficient  C^ 
Tesseral  harmonic  coefficient  S^ 

S33 

1 .4343E-6 

C44 

-0. I846E-6 

Tesseral  harmonic  coefficient  C. , 

44 

S44 

0.2508E-6 

Tesseral  harmonic  coefficient  S. . 

44 

C31 

2.0192E-6 

Tesseral  harmonic  coefficient 

S3! 

0.2278E-6 

Tesseral  harmonic  coefficient  S^j 

C42 

0. 3444E-6 

Tesseral  harmonic  coefficient  C ^ 

S42 

0.7021E-6 

Tesseral  harmonic  coefficient  S^ 

ERAD 

6378.163  km 

Earth’s  mean  equatorial  radius 

E OMEGA 

7.2921 15I47E-5  rad  sH 

Earth's  rotation  rate 

PREC 

6.079E-12  rad  s_l 

. . 3 -2 

Earth's  precession  rate 

SUNMU 

1 . 327 1 27E+1 1 km s 

Sun's  gravitational  constant 

SELMU 

4902.756  km3  s"2 

Moon's  gravitational  constant 

Data  in  /PERTURB/  - 

Variable  name 

Value 

Explanation 

HALF CD 

1.1 

Product  JCp  where  Cp  is  the 
drag  coefficient. 

RATE 

1.1 

Ratio  of  the  angular  velocity  of 
the  atmosphere  to  that  of  the 
earth. 

I 
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SUBROUTINE  COTOEL 


Summary 


Language 


Author 


- The  subroutine  derives  the  standard  elliptic  (osculat- 
ing) orbital  elements  of  an  earth  satellite,  given  its 
position  and  velocity  components. 

- ASA  Fortran  (Standard  Fortran  4) . 

- R.H.  Gooding  (May  1976). 


Subroutine  statement  - SUBROUTINE  COTOEL  (X,  Y,  Z,  XDOT,  YDOT,  ZDOT,  EMU,  A, 

E,  ORBINC,  RANODE,  ARGPER , EM,  EN) . 


Input  arguments 


X,  Y,  Z 


XDOT,  YDOT, 


Output  arguments 


ORBINC 


RANODE 


ARGPER 


Use  of  COMMON 


Source  deck 


Local  storage  used 


Geocentric  cartesian  components,  (x,y,z),  of  the  satellite, 
x towards  the  vernal  equinox  and  z towards  the  north 
pole. 

Velocity  of  the  satellite,  (x,y,z),  measured  in  the 
same  coordinate  system. 

Earth's  gravitational  constant,  u . 


Semi -major  axis,  a , in  the  same  units  as  x,  y,  z 
Eccentricity,  e . 

Orbital  inclination,  i . 

Right  ascension  of  the  ascending  node,  0 . 

Argument  of  perigee,  a>  . 

Mean  anomaly,  M . 

Mean  motion,  n . 

- None. 

- 30  cards,  including  four  comment  cards  (ICL  code). 

- 13  real  variables. 


Subordinate  subprograms  - The  subroutine  TRINV. 


Explanation  - Although  the  subroutine  essentially  transforms  from  the 

six  components  of  the  position  and  velocity  of  a satellite  to  its  six  orbital 
elements,  a seventh  input  argument  permits  an  arbitrary  value  of  the  constant  y 


o\ 

r*. 
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Co  be  used;  Che  sevench  output  argument,  n , is  derived  from  a and  p by  Che 
2 3 

relation  n a ■ p (Kepler's  third  law).  This  means  that  the  subroutine  may  be 
used  very  generally;  e.g.  for  a planet,  taking  the  sun's  p and  interpreting  A 
as  celestial  longitude.  Units  of  time  and  distance  are  arbitrary  but  must,  of 
course,  be  consistent;  all  angles  are  in  radians. 

The  subroutine  has  been  written  to  give  maximum  accuracy.  All  angles  are 
. derived  from  knowledge  of  both  sine  and  cosine,  and  in  such  an  order  that  there 

is  no  difficulty  near  the  singularities  at  e - 0,  i - 0 and  i = n . The 
general  method  is  that  of  Brouwer  and  Clemence  (section  27  of  chapter  1). 

The  first  quantities  to  be  derived  are  a,  e and  the  eccentric  anomaly, 

E . (NB  The  Fortran  variable  E refers  to  the  eccentricity  and  not  the  eccentric 
anomaly.)  These  come  from  the  relations 


a 


2 V 


e cos  E - 


rV 


- 1 


and 


e sin  E 


(xx  + yy  + zz) 

. iJ.  . 

(ua)* 


. , 2 ^ 2 2*J  , „2  .2  .2  .2 

where  r ■ (x  ♦ y ♦ z ) and  V ■*  x + y + z 

Then  M follows  at  once,  since 


M ■ E - e sin  E 


Note  that  if  e is  close  to  zero  E is  ill-determined,  reflecting  the 
indeterminancy  of  perigee  - and  in  fact  if  e - 0,  E is  set  to  0 . This  does 
not  matter  at  all,  since  whatever  position  of  perigee  is  specified  by  the  value 
taken  for  E , the  value  of  M is  fully  consistent  with  it. 

The  values  of  i,  0 and  u are  derived  from  the  basic  matrix  relation 

P P P \ 
x y z\ 

Qx  Qy  Qz  I 

Rx  Ry  \l 


f COM  U 

sin  u 

°\ 

0 

0 

\ / cos  A 

sin  A 

°\  / 

-sin  a) 

COS  It) 

° 

0 

cos  i 

8 in 

i]  [-sin  A 

cos  A 

0 - 

» 

0 

>/ 

lo 

-sin  i 

cos 

J \ 0 

0 

1/  \ 

mm 
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where  - (x/r)  cos  E - x(a/p) ^ sin  E , 

Qx  - 0 “ e2)  ^ [(x/r)  sin  E + x(a/p)*(cos  E - e)]  , 

Rx  - (ya(l  - e2))  ^ (yz  - zy ) , 

and  similarly  for  P , P , Qv,  Q , R and  R . 

y z v z y z 

Thus  i and  ft  are  derived  from 


and 


with  an  indeterminacy  in 
sin  i - 0). 


sin  i sin  ft  - R , 
x * 

- sin  i cos  ft  ■ R , 

y 

cos  i • R , 
z 

ft  when  i is  close  to  0 or 


IT 


(ft  is  set  to  0 if 


Finally  u is  determined  from 


and 


sin  i cos  u ■ Q 

z 

sin  i sin  u = P 

z 


where  we  have  to  be  sure  that  no  difficulty  arises  over  the  e-singularity  or 
either  of  the  i-singularities . 

There  is  no  trouble  near  the  e-singularity  since  the  quantities  e cos  E 

and  e sin  E are  used  (through  Qz  and  P^),  in  the  derivation  of  w , in  the 

same  ratio  as  in  the  derivation  of  E , so  that  E + to  is  always  correct. 

(When  E is  set  conventionally  to  0 because  e - 0 , e cos  E is  set  to  1 (I) 

to  ensure  the  correct  ratio  between  Q and  P .) 

z z 

In  the  same  way,  w has  to  be  compatible  with  ft  near  the  i-singularity , 

and  this  is  automatic  when  z and  z are  not  both  exactly  zero,  since  ft  and 

u both  depend  on  the  ratio  of  these  two  quantities.  When  z and  z are  both 

exactly  zero,  the  correct  value  of  u>  is  achieved  by  replacing  P and  Q by 

z z 

Py  and  Qy  , with  an  additional  factor  to  ensure  that  u is  in  the  correct 
quadrant. 
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SUBROUTINE  DR 


Stannary  - The  subroutine  computes  the  direction  cosines  of  the 

nominal  ABM  firing  direction  required  to  obtain  a 
specified  drift  rate  after  ABM  burn. 


Language 

Author 

Subroutine  statement 
Input  arguments 


- 1900  Fortran. 

- M.D.  Palmer  (May  1975). 

- SUBROUTINE  DR  (AD) 


AD  The  drift  orbit  semi  major  axis  corresponding  to  the 

required  drift  rate  (km). 

Use  of  COMMON  - Certain  arguments  in  the  common  blocks  /ABM/,  /CNST/ 

and  /CONST/  are  used  as  follows 


Input  arguments  in  /ABM/  - 

The  satellite's  transfer  orbit  geocentric  position, 

-b  " ^x»  y’  z^’  and  vel°city,  vfe  - (x,  y,  z),  compon- 
ents at  the  ABM  burn  point  (km,  km  s *)• 

The  satellite's  radial  distance,  r , and  velocity, 

i 

, at  the  ABM  burn  point  (km,  km  s ) . 

Direction  cosines  of  the  unit  vector,  h , normal  to 
the  required  drift  orbit  plane. 

Nominal  ABM  velocity  increment  (km  s *). 

Input  argument  in  /CNST/  - 

PI  ir  . 

PI2  ir/2  . 

Input  argument  in  /CONST/  - 

3 -2 

EMU  Earth's  gravitational  constant  (km  s ). 

Output  arguments  in  /ABM/  - 

SX,  SY,  SZ  Direction  cosines  of  the  required  ABM  firing  direction. 

Source  deck  - 56  cards.  (ICL  code). 


PT(3)  I 
VT  (3)  J 

R,  VB 

HX,  HY,  HZ 

DV 


Local  storage  used  - 37  real  variables. 


Subordinate  subprograms  - None. 
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Explanation  - The  subroutine  calculates  the  direction  cosines  of  the 

ABM  firing  direction  required  to  produce  a specified  drift  rate  after  ABM  burn. 
This  parameter  is  passed  to  the  subroutine  as  the  equivalent  drift  orbit  semi 
major  axis. 

The  drift  orbit  velocity  required  after  ABM  burn  is  calculated  together 
with  the  angle  (0)  between  the  transfer  and  drift  orbit  velocity  vectors  and  the 
angle  ($)  between  the  transfer  orbit  velocity  vector  and  the  unit  vector  h . 

If  the  ABM  velocity  increment  is  inadequate,  ie  <J>  + 0 < 90°  , a STOP  ABM  DV 
instruction  is  executed. 


Let  the  direction  cosines  of  the  drift  orbit  velocity  vector,  the  unit 

a 

vector  h and  the  transfer  orbit  velocity  vector  be  (A^,  Ay,  1^),  (h^,  hy,  hz> 

and  (b  , b , b ) respectively.  Then 
x y z 

b 1 + b A + b A - cos  0 , 

x x y y z z ’ 


and 


2 2 2 

A + A + A 
x y z 


1 . 


h 1 + h A + h i * 0. 

xx  y y z z 


Therefore 

A - ■ 

z 

and 

b A + b A 
xx  y y 

ie 

k - d1 

or 

AA. 

Setting 

A 

X 

we  have 

A2  + A2  + [h2A 
x y L x : 

ie 

L x x y yJ  z 

b [h  A + h l 1 
zL  x x y yJ 

h 

z 

k * k - t*k 


cos  0 


cos  0 


A^  ■ [cos  0 - BAy]/A  , 


,2.2-.  ,.2 


fcos  8 _ B T 

LA  A VI 


1 + -4 


h2 

_x 

2 


2 f + A 

L A A yJ  y y 


1 + Hr 

hf 


1 
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or 

af.^  + bi.  ♦ c ■ 0 . 

y y 

2 2 
If  b - 4ac  is  negative,  the  program  obeys  a STOP  90  instruction.  If  b > 4ac 

the  equation  has  two  possible  solutions.  These  are  evaluated  and  the  correspond- 
ing values  of  l and  determined.  The  drift  orbit  eccentricity  for  each 

solution  is  calculated  and  the  direction  cosines  of  the  ABM  firing  direction 
determined  for  the  solution  which  gives  the  minimum  drift  orbit  eccentricity. 
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Summary 


SUBROUTINE  DVST 

- The  subroutine  calculates  the  delta-velocity  required 
fo_-  a geosynchronous  satellite  to  acquire  station 
after  ABM  burn.  Both  easterly  and  westerly  drifts 
are  examined  and  the  one  giving  the  smaller  total 
delta-velocity  selected. 


Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (January  1974). 


Subroutine  statement  - SUBROUTINE  DVST (NO) 
Input  argument 


NO  The  case  number.  If  NO  ■ 0,  no  line  printer  output 

is  produced. 

Use  of  COMMON  - Certain  quantities  in  the  common  blocks  /CNST/, 

/STATION/  and  /STATS/  are  used  as  follows s- 

Input  arguments  in  /CNST/  - 

ENO  Earth's  rotation  rate  (rad/d) 

RAD  Conversion  factor,  degrees  to  radians. 

Input  arguments  in  /STATION/  - 

X,  Y,  Z Satellite's  geocentric  position  components, 

r^  * (x,  y,  z) , at  ABM  burn  (km). 

XDOT,  YDOT,  ZDOT  Satellite's  geocentric  velocity  components, 

v,  * (x,  y,  z) , immediately  after  ABM  burn  (km  s *). 
-a 


MJDB 

TB 


TR 

Cl 


PHI 


Modified  Julian  day  number  corresponding  to  the  ABM 
burn  time. 

Time  of  ABM  burn,  in  fractions  of  a day,  relative 
to  MJDB. 

The  time  allowed  for  tracking  and  orbit  determination 
after  ABM  burn  (days). 

Cl  + TR  is  the  total  time  available  for  station 
acquisition  (days). 

Station  longitude  (rads  E of  Greenwich). 
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XIL,  XIH 


83 


Minimum  and  maximum  permitted  values  of  drift  orbit 
inclination  (rad) . 


Input  arguments  in  /STATS/  - 


RAS,  DEC 

Right  ascension  and  declination  of  the  actual  ABM  fir- 
ing direction  (rad) . 

Output  arguments  in  /STATS/  - 

DV 

Total  delta-velocity  required  for  station  acquisition 
(ms  1 ) . 

XINCD 

Drift  orbit  inclination  (degree). 

BO 

Right  ascension  of  the  ascending  node  of  the  drift 

orbit  (degree) . 

XLOND 

Satellite's  longitude  at  the  ABM  burn  point  (degree  E) 

A 

Drift  orbit  semi  major  axis  (km). 

E 

Drift  orbit  eccentricity. 

Source  deck 

- 49  cards  (ICL  code). 

Local  storage  used 

- 4 real  array  elements;  20  real  variables  and  1 integer 

variable. 

Subordinate  subprograms 

- The  subroutines  COTOEL  and  TRINV  and  the  functions 

ANGLE,  EAFKEP  and  SIDANG. 

Explanation  - The  subroutine  uses  linearised  equations  to  calculate 

the  total  delta-velocity  required  for  a synchronous  satellite  to  acquire  station. 
Both  easterly  and  westerly  drifts  are  examined  and  the  one  requiring  the  smaller 
delta-velocity  selected. 

The  value  of  the  earth's  gravitational  constant  is  taken  to  be 
3 -2 

398616.82  km  s which  takes  into  account  the  effects  of  the  earth's  second 
zonal  harmonic  • 

COTOEL  is  called  to  find  the  osculating  elements  of  the  drift  orbit 
immediately  after  ABM  burn.  The  satellite  longitude  at  the  ABM  burn  point  is 
calculated.  If  NO  0 , the  case  number,  the  ABM  firing  attitude,  drift  orbit 
elements  and  the  longitude  of  the  ABM  burn  point  are  output  to  the  lineprinter. 
The  inclination  is  compared  with  the  maximum  and  minimum  permitted  values  and, 
if  necessary,  the  delta-velocity  required  to  place  the  inclination  within  these 
bounds  is  determined. 
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Both  easterly  and  westerly  drifts  are  now  examined.  In  each  case,  the 
necessary  drift  rate  is  computed.  If  the  actual  drift  rate  is  inadequate  or  in 
the  wrong  direction  a manoeuvre  is  performed.  It  is  assumed  that  this  and  the 
final  stopping  manoeuvre  are  also  used  to  reduce  the  orbital  eccentricity. 
Finally  the  total  velocity  increment  needed  to  remove  any  residual  eccentricity 
is  calculated.  If  NO  ^ 0 , one  line  is  printed  giving,  for  each  direction,  the 
drift  rate  and  the  velocity  increments  required. 
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Summary 


Language 

Author 

Function  statement 
Input  arguments 
EM 
ECC 

Output  function 
EAFKEP 

Use  of  COMMON 

Source  deck 

Local  storage  used 

Subordinate  subprograms 

Explanation 

Improved  approximations 
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FUNCTION  EAFKEP 

- The  function  solves  Kepler's  equation;  i.e.  it  provides 
the  eccentric  anomaly  of  a celestial  body,  E , given 
the  orbital  eccentricity,  e , and  mean  anomaly,  M . 
Kepler's  equation  is  M - E - e sin  E . 

- ASA  Fortran  (Standard  Fortran  4). 

- R.H.  Gooding  (May  1976). 

- FUNCTION  EAFKEP  (EM,  ECC). 

Mean  anomaly,  M (radians). 

Eccentricity,  e . 

Eccentric  anomaly,  E . 

- None. 

- 13  cards,  including  3 comment  cards  (ICL  code). 

- 3 real  variables. 

- None. 

- A first  approximation  to  E is  given  by  E|  - M . 
are  given  by  Newton's  method;  thus 


The  process  ends  after  three  iterations  if  e < 0.003  , and  otherwise  after  five 
This  ensures  sufficient  accuracy  at  all  times,  while  making  the  subroutine 
independent  of  the  word-length  of  the  computer.  (If  the  magnitude  of  |e^  + j “ E^ 
were  used  as  a criterion  for  convergence,  the  numerical  value  it  was  compared 
with  would  have  to  vary  from  computer  to  computer.) 
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Summary 


Language 

Author 

Subroutine  statement 
Input  argument 
GAM 

Use  of  COMMON 


SUBROUTINE  FPA 

- The  subroutine  finds  the  direction  cosines  of  the 
nominal  apogee  boost  motor  (ABM)  firing  direction 
required  to  obtain  a specified  flight  path  angle  at 
ABM  burn. 

- 1900  Fortran. 

- M.D.  Palmer  (May  1975). 

- SUBROUTINE  FPA  (GAM) 


The  required  flight  path  angle,  y (rad). 

- Certain  quantities  in  common  block  /ABM/  are  used 
as  follows 


Input  arguments  in  /ABM/  - 


PT(3) 

VT(3) 


RB,  VB 

HX,  HY,  HZ 

DELTAV 

Output  arguments 
SX,  SY,  SZ 

Source  deck 


The  satellite's  transfer  orbit  geocentric  position, 

r - (x,  y,  z),  and  velocity,  v - (x,  y,  z) , 
b . ~ “ - 1 

components  at  the  ABM  burn  point  (km,  km  s ) . 

The  satellite's  radial  distance  and  velocity  at  the 
ABM  burn  point  (km,  km  a *). 

Direction  cosines  of  the  unit  vector  (h)  normal  to 
the  required  drift  orbit  plane. 

Nominal  ABM  velocity  increment  (km  s * ) . 


Direction  cosines  of  the  required  ABM  firing 
direction. 

- 24  cards  (ICL  code) . 

- 12  real  variables. 


Local  storage  used 
Subordinate  subprograms  - None. 

A 

Explanation  - Let  r^  be  a unit  vector  along  the  radius  vector  at 

the  ABM  burn  point  and  define  the  unit  vector  u such  that  u,  h,  r^  form  a 
right  handed  orthogonal  set  (see  Fig  5).  Then 

A A A 

u - I*  * * . 


r \ 


179  Appendix  D 87 

The  flighc  path  angle  (y)  is  measured  from  u in  the  plane  defined  by 

A A A 

u and  rfe  . If  vd  is  a unit  vector  along  the  drift  orbit  velocity  vector 
vd  , then 


and 


■ u cos  y + r^  sin  y 

|vd|  “ |Yb  . vd|  + [6v2  - (v2  - |vb  . vd | 2)]i  , 


where  6v  is  the  nominal  ABM  velocity  increment.  If  the  argument  of  the 
square  root  is  negative,  a STOP  FPA  instruction  is  executed.  The  direction 
cosines  of  the  nominal  ABM  firing  direction,  (s^,  sy,  s ),  are  found  from  the 
equation 


s * 

X 

[M 

1*.  - 
X 

iv 

M 

1 6v  , 

sy  " 

L|vd' 

In  - 
1 y 

K' 

j 6v 

S * 

z 

h' 

k - 

z 

W 

^6v 

and 


where  (£,£,£)  are  the  direction  cosines  of  the  drift  orbit  velocity  vector 
x y z 

and  (b  , b , b ) are  the  direction  cosines  of  the  transfer  orbit  velocity  vector, 
x y z 
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Summary 


Language 


Author 


SUBROUTINE  INTPTB 

- The  subroutine  calculates  the  cartesian  components  of 
the  geocentric  position  and  velocity  of  an  earth  satel- 
lite at  a given  time  by  interpolation  between  two 
points  in  the  orbit.  These  points  are  defined  by 
position,  velocity  and  acceleration  components  at 
specified  times. 

- 1900  Fortran. 

- M.D.  Palmer  (July  1974) 


Subroutine  statement  - SUBROUTINE  INTPTB  (TINT,  TP,  TIMET,  PPVA,  PT,  VT) . 


Input  arguments 


TIMET 


PPVA (9) 


Output  arguments 


PT(3) 


VT(3) 


Use  of  COlfclON 


Input  arguments 


PP(3) 


XACCT(3) 


Time,  t , at  which  the  interpolation  is  required. 

Time,  tj  , of  the  point  prior  to  the  interpolation. 

Time,  t£  , of  the  point  after  the  interpolation. 

Satellite's  geocentric  position,  velocity  and  accelera- 
tion components  at  time  tj  (Xj,yj,*j  etc.). 

Satellite's  geocentric  position  components  at  time  t 

(x.y.z). 

Satellite's  geocentric  velocity  components  at  time  t 
(x,y,z). 

- Certain  arguments  in  the  common  block  /CINTEG/  are  used 
as  follows: 


Satellite's  geocentric  position  components  at  time  t2 

(*2 »y  2»  z2^  * 

Satellite's  geocentric  velocity  components  at  time  t2 

(Xj.yj.^* 

Satellite's  geocentric  acceleration  components  at  time 

t2  (x2,y2,l2). 
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Local  8 tor age  used  - 23  real  variables  and  I integer  variable. 

Subordinate  subprograms  - None. 

Source  deck  - 40  cards  (ICL  code). 

Explanation  - The  subroutine  obtains  a satellite's  geocentric 

position  and  velocity  components  at  a given  time  by  interpolation  between  two 
points  in  the  orbit.  These  will  normally  be  adjacent  steps  in  a numerical 
integration. 

These  equations  used  for  interpolation  are: 
f(t)  - q3[o  ♦ 3p  ♦ 6p2)f  (t])  + hp(l  + 3p)f'(t1)  + p2?'^,)! 

+ 3q  + 6q2)f(t2)  — hq ( 1 + 3q)f  ' (t2) 

for  position  components  and 

f'(t)  - 30h_,p2q2[f(t2)  - f(tj)]  + q2(l  + 5p) ( 1 - SpJf’Ct,) 

+ p2(l  + 5q)(l  - 3q)f '(t2) 

+ | Pq2(2  - 5p)f " (t j ) - ^ p2q(2  - 5q)f"(t2) 

for  velocity  components,  where  h ■ t2  - t j , p ■ (t  - tj)/h,  q ■ I - p and 
f(t^),  f'(t.)  and  f"(t^)  are  the  position,  velocity  and  acceleration  components 
at  time  t^  . 
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Summary 

Language 

Author 


SUBROUTINE  MATADD 

- The  subroutine  adds  two  matrices  of  order  m x n . 

- 1900  Fortran. 

- M.D.  Palmer  (July  1974) 


Subroutine  statement  - SUBROUTINE  MATADD  (A,  B,  C,  M,  N) . 


Input  arguments 


Output  arguments 


Matrices  of  order  M x N 


Number  of  rows. 


Number  of  columns. 


Use  of  COMMON 
Source  deck 
Local  storage  used 


Matrix  of  order  M x N (C  - A + B) . 


- 8 cards  (ICL  code). 


- 2 integer  variables. 


Subordinate  subprograms  - None. 


Explanation 

C - A + B.  Each  element 
elements  of  A and  B . 


■ The  subroutine  performs  the  matrix  operation 
of  C is  obtained  by  the  addition  of  the  corresponding 
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SUBROUTINE  MATMUL 

- The  subroutine  performs  matrix  multiplications. 

- USA  Standard  Fortran  (USAS  X3.9  - 1966). 

- G.E.  Cook  (May  1970). 

- SUBROUTINE  MATMUL  (EMI,  EM2,  EM3,  II,  KK,  JJ,  IA,  IB). 


EMI 

EM2 

II 

KK 

JJ 

IA 

IB 


Output  arguments 


EM3 


Matrix  of  dimension  (II, KK). 

Matrix  of  dimension  (KK,JJ). 

Actual  number  of  rows  in  EMI  . 

Actual  number  of  columns  in  EMI  and  rows  in  EM2  . 

Actual  number  of  columns  in  EM2  . 

Maximum  number  of  rows  in  EMI  as  given  in  the 
dimension  statement  of  the  calling  segment. 

Maximum  number  of  rows  dimensioned  for  EM2  in  the 
calling  segment. 

Matrix  of  dimension  (II, JJ). 


Use  of  COMMON  - None. 

Source  deck  - 13  cards,  including  4 conment  cards  (ICL  code). 

Local  storage  used  - 3 integer  variables. 


Subordinate  subprograms  - None. 
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SUBROUTINE  OFPA 

Summary  - The  subroutine  optimises  the  drift  orbit  flight  path 

angle  at  ABM  burn  such  that  the  total  delta-velocity 
required  for  station  acquisition  and  circularisation 
of  the  drift  orbit  is  a minimum. 


Language  - 1900  Fortran. 

Author  - M.D.  Palmer  (September  1975). 

Subroutine  statement  - SUBROUTINE  OFPA  (STAC). 


Input  argument 

STAC  The  name  of  the  subroutine  to  be  used  to  compute  the 

delta-velocity  required  for  station  acquisition. 

Use  of  COMMON  - Certain  quantities  in  conmon  blocks  /ABM/,  /CNST/, 

/STATION/  and  /STATS/  are  used  as  follows :- 

Input  arguments  in  /ABM/  - 


PT(3)1 

VT(3)J 


The  satellite's  geocentric  transfer  orbit  position, 
Ib  - (x,  y,  z),  and  velocity,  vfa  - (x,  y,  z) , 
components  at  the  ABM  burn  point  (km,  km  s 1 ). 


RB,  VB 


HX,  HY,  HZ 

DV 

DRA 

SDIN 


The  satellite's  radial  distance,  r.  , and  the 

D 

transfer  orbit  velocity,  v,  , at  the  ABM  burn  point 

-1  ® 

(km,  km  s ). 

Direction  cosines  of  the  unit  vector,  h , normal  to 
the  required  drift  orbit  plane. 

The  nominal  ABM  velocity  increment  (km  s *) 

Right  ascension  of  the  ascending  node  of  the  required 
drift  orbit  (rad). 

Sine  i,  , where  i,  is  the  desired  drift  orbit 
d d 

inclination. 


Input  argument  in  /CNST/  - 


PI 


IT  . 


RAD 


Conversion  factor,  degrees  to  radians 
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Input  argument  in  /STATS/  - 


Total  velocity  increment  required  for  station 
acquisition  (ms  ) . 


Output  arguments  in  /ABM/  - 
SI,  S2,  S3  I 


Direction  cosines  of  the  required  ABM  firing 
direction. 


Output  arguments  in  /STATION/  - 


X,  Y,  Z 


Source  deck 


Local  storage  used 


Satellites  geocentric  position  components  at  the 
ABM  burn  time  (km) . 

Satellite's  geocentric  drift  orbit  velocity  com- 
ponents, v^  - (x,  y,  z),  immediately  after  ABM 
burn  (km  s~* ) . 

- 113  cards  (ICL  code). 

- 37  real  variables. 


Subordinate  subprograms  - The  subroutines  ANST,  COTOEL,  DVST,  and  TRINV  and 

the  functions  ANGLE,  EAFKEP  and  SIDANG. 

Explanation  - The  subroutine  optimises  the  drift  orbit  flight  path 

such  that  the  total  velocity  increment  required  for  station  acquisition  is  a 

minimum.  The  value  of  the  earth's  gravitational  constant  is  taken  to  be 
3 -2 

398616.82  km  s which  takes  into  account  the  effects  of  the  earth's  second 
zonal  harmonic,  . 

The  transfer  orbit  inclination  (i^)  and  the  required  nodal  rotation  (AD) 
are  calculated.  The  angle,  Aa  , between  the  transfer  and  drift  orbit  planes  is 
then  obtained  using  the  equation 

cos  Aa  ■ sin  i,  sin  i cos  AD  + cos  i,  cos  i , 
at  d t 

where  i^  is  the  required  drift  orbit  inclination.  The  angle  <J>t  between  the 
transfer  orbit  radius  and  velocity  vectors  is  given  by 

cos  *t  - Eb  • Yb  , 

where  r^  and  v^  are  unit  vectors  along  the  transfer  orbit  radius  and 
velocity  vectors  respectively. 
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The  angle  between  the  transfer  and  drift  orbit  velocity  vectors  is 


given  by 


cos  - cos  <j>t  cos  + sin  <J>t  sin  <J>d  cos  A a , 


where  <j>  is  the  angle  between  the  drift  orbit  radius  and  velocity  vectors, 
d 

If  the  nominal  ABM  velocity  increment  is  6v  , 

6v2  - vd  + Vb  ‘ 2Vb  C0S  * ■ 

and  so  given  iji  , v^  can  be  determined.  The  subroutine  uses  an  iterative 

procedure  to  find  the  direction  cosines  (Jl  , 1 , n ) of  the  drift  orbit 

x y z 

velocity  vector.  Let  (b  , b , b ) and  (h  , h , h ) be  the  direction  cosines 

x y z x y z 

of  the  transfer  orbit  radius  vector  at  the  ABM  burn  point  and  the  normal  to  the 
drift  orbit  plane. 


h.v,  - h n + h l + h t ■ 0 , 

-d  xx  yy  zz 


2 2 2 

n + n + n = i , 

x y z 


r,  .v,  ■ bH  + b l + b l 

~b  -d  xx  y y z z 


cos  $ . 


I * [-  h l + h 1 1/h 
x 1 y y z zJ'  x 


cos  ■ bJl  + b n - b [h  n + h n 1/h  , 

d y y z z xL  y y z z x 


from  which  we  obtain 


An  + Bn  ■ cos  , 
y z d 


It  follows  that 


" Ccos  ♦d  “ BAJ/A 


0 0 h2  2h  h n n 

L i}  + o2  _z  + , + y z y z 

2 z y .2  1 2 

h n 

xj  L x J x 


j 
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and  substituting  for  i 


hht 
y z z 


cos  4> . 


-k] 


These  equations  are  solved  to  obtain  two  values  of  Hz  and  corresponding 

values  of  l and  £,  . The  product  v,  . v,  is  formed  for  each  solution  and 
y x -d  -b 

the  components  of  v,  are  computed  for  the  scalar  product  which  is  closest  in 
— d 

value  to  cos  ij>  . A subroutine  is  then  called  to  compute  the  total  velocity 
increment  required  for  station  acquisition. 

In  practice,  is  set  initially  to  80  degrees  and  the  delta-velocity 

computed  and  stored.  <J>d  is  then  incremented  by  6 , (4  degrees  initially)  and 
the  procedure  repeated  until  the  delta-velocity  passes  through  a minimum. 

6a^  is  then  divided  by  four  and  the  process  is  repeated  until  a minimum  is 
found  with  6<f>d  < 0.001  radian.  The  appropriate  set  of  direction  cosines  are 
found  as  above  and  the  direction  cosines  of  the  nominal  ABM  firing  direction 
are  obtained  using  the  equation 

6vSl  - |vd|Hx  - v^  , 

6vs2  - lYdl*y-vby  , 


53  “ " Vb2  * 
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FUNCTION  SIDANG 


Summary 

Language 

Author 

Function  statement 
Input  arguments 
MJD 
TIME 

Output  function 


- The  function  calculates  the  modified  sidereal  angle 
corresponding  to  a given  date/time  (t). 

- A.S.A.  Fortran  (Standard  Fortran  4). 

- R.J.  Tayler  (September  1969) 

- FUNCTION  SIDANG  (MJD,  TIME) . 

Modified  Julian  day  number 

Time  (fraction  of  a day)  such  that  t is  given  by 
MJD  + TIME. 


SIDANG 

Use  of  COMMON 
Source  deck 
Local  storage  used 
Subordinate  subprograms 


The  modified  sidereal  angle  (radians) 

- None. 

- 11  cards,  including  5 comment  cards  (ICL  code) 

- None. 

- None. 


Explanation  ~ The  modified  sidereal  angle  is  a quantity  which 

differs  from  the  more  familiar  sidereal  time  only  because  of  the  choice  of  a 
non-standard  reference  direction.  This  reference  direction  lies  in  the  plane 
of  the  true  equator  and  is  given  by  rotation  of  the  true  equinox,  by  an  amount 
equal  to  the  precession  and  nutation  in  right  ascension  since  1950.0  but  in 
the  opposite  sense. 

In  degrees , the  formula  for  modified  sidereal  angle  is 

0 = 100°. 075542  + 360°. 9856 12288  (d  - 33282.0) 

ie  the  origin  has  been  adjusted  slightly  from  1950.0,  which  occurred  at 

MJD  33281.9234.  However,  SIDANG  operates  by  calculating  tha  angle  in  revolu- 
# I 

tions,  as 

. 

0.277987616  + 0.00273781191  (MJD  - 33282)  + 1.00273781191  TIME 

with  an  integral  number  of  revolutions  already  dropped,  and  then  converting  to 
radians  (after  dropping  any  remaining  full  revolutions)  by  multiplying  by 
6.2831853072. 
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SUBROUTINE  STATIS 

Summary 

“ The  subroutine  produces  statistical  information  about 

the  orbital  elements  of  the  drift  orbit  and  the  delta- 
velocity  required  for  station  acquisition. 

Language 

- 1900  Fortran. 

Author 

- M.D.  Palmer  (April  1976). 

Subroutine  statement 

- SUBROUTINE  STATIS  (RNC) 

Input  argument 

- 

RNC 

On  the  first  n calls,  n being  the  number  of 
cases,  RNC  *>  0 . On  the  n + 1 th  call,  RNC  * 1/n. 

(see  explanation). 

Use  of  COMMON 

Certain  quantities  in  common  block  /STATS/  are  used 

as  follows:- 

Input  arguments 

in  /STATS/  - 

RA  1 
DEC  J 

1 

Right  ascension  and  declination  of  the  ABM  firing 
direction  (degree). 

DV 

Delta-velocity  required  for  station  acquisition 
(m/s) . 

01 

Drift  orbit  inclination  (degree). 

BO 

Right  ascension  of  the  ascending  node  of  the  drift 
orbit  (degree). 

XLON 

Satellite  longitude  at  ABM  burn  (degrees  east  of 

Greenwich) . 

A 

Drift  orbit  semi  major  axis  (km). 

E 

Drift  orbit  eccentricity. 

THEB 

Solar  aspect  angle  at  ABM  burn  (degree) . 

ELEV 

Spin  axis  elevation  angle  at  ABM  burn  (degree). 

NTHEB 

The  number  of  cases  for  which  the  solar  aspect  angle 

constraint  has  been  violated. 

Source  deck 

- 147  cards  (ICL  code). 

Local  storage  used  - 22  integer  variables  and  41  real  variables. 
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Subordinate  subprograms  - None . 

Explanation  - On  the  ith  call  (i  ■ l,...,n)  STATIS  performs  the 

following  operations  on  each  quantity  in  the  common  block  /STATS/, 


X.  “ X.  + x. 
1 l-l  i 


Y.  ■ Y.  . ♦ xt 

l l-l  i 


where  Xq  ■ Y^  ■ 0 . 

The  satellite  longitude  is  checked  to  ensure  that  it  is  in  the  correct 
quadrant,  and  the  delta-velocity  required  for  station  acquisition  is  stored  in 
a form  suitable  for  subsequent  output  as  a histogram. 


When  the  n simulations  are  complete,  STATIS  is  called  with  RNC  - 1/n 
and  the  mean  and  standard  deviations  of  the  above  quantities  are  calculated 
using  the  formulae :- 


and 


This  data  is  output  on  the  linepr inter,  together  with  a histogram  of 
delta-velocities  and  a comment  specifying  the  number  of  cases  for  which  the 
solar  aspect  angle  constraint  has  been  violated. 
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SUBROUTINE  XROT 


Summary 


Language 

Author 

Subroutine  statement 
Input  arguments 
THETA 


Y 

Z 


- Rotation  of  a rectangular  coordinate  system  through 
a specified  angle  in  a single  plane. 

- 1900  Fortran. 

- Diana  W.  Scott  (April  1969). 

- SUBROUTINE  XROT  (THETA,  Y,  Z) . 


Angle  (in  radians)  through  which  system  is  to  be 
rotated  about  the  x axis. 

Cartesian  y coordinate. 

Cartesian  z coordinate. 


Output  arguments 


Y 

Z 

Use  of  COMMON 

Source  deck 

Local  storage  used 

Subordinate  subprograms 

Explanation 


Rotated  y coordinate. 

Rotated  z coordinate. 

- None. 

- II  cards,  including  2 comment  cards  (ICL  code). 

- 3 real  variables. 

- None. 

- y and  z are  replaced  by  (y  cos  6 + z sin  6)  and 


(z  cos  0 - y sin  0)  respectively.  The  subroutine  can  of  course  be  used  for 
rotations  about  the  y or  z axes  as  well,  eg  to  rotate  about  the  y axis,  write 
CALL  XROT  (THETA,  Z,  X). 
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